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Abstract 


A  radioactive  fallout  code,  HYDRA,  has  been  developed 
chat  uses  real  winds  to  determine  dose  rate  contours  on  the 
ground.  These  winds  are  calculated  from  spectral  coeffi¬ 
cients  derived  by  the  National  Meteorological  Center. 

HYDRA  models  the  radioactive  dust  cloud  as  a  set  of  pancake 
clouds,  each  cloud  representing  a  different  particle  size 
group.  Each  particle  size  group  is  transported  to  the 
ground  through  a  discretely  layered  atmosphere  using 
McDonald-Da vies  fall  mechanics.  Dose  rate  contours  are 
determined  by  smearing  the  cloud  activity  along  the  ground 
as  the  cloud  descends. 

HYDRA  was  evaluated  against  another  variable  wind 
fallout  code,  REDRAM.  HYDRA's  multiple  pancake  cloud 
model,  as  compared  to  the  single  pancake  cloud  model  used 
in  REDRAM,  produced  the  same  contour  patterns  for  yields 
between  10  and  2000  kilotons.  Outside  this  range,  however, 
dose  rates  along  the  hotline  were  lower  and  the  contours 
wider  for  HYDRA's  multiple  pancake  cloud  model. 

The  Colarco  coefficients  used  in  REDRAM  were  found  to 
consistently  underestimate  particle  size.  Variations  were 
as  high  as  six  percent  for  yields  ranging  from  100  to  1000 
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kilotons.  HYDRA  did  not  use  Colarco's  approximation. 

HYDRA  successfully  produced  dose  rate  contours  at  nine 


global  loca tions\ 


A  COMPUTER  FALLOUT  MODEL  USING  VARIABLE 
WINDS  FOR  OPERATIONAL  TYPE  STUDIES 


I .  Introduction 

Radioactive  Cloud  Formation  and  Fallout 

When  a  nuclear  weapon  is  detonated,  a  tremendous 
amount  of  energy  is  released,  mostly  as  X-ray  photons.  The 
weapon  material  and  casing  are  initially  subjected  to  tem¬ 
peratures  as  high  as  ten  million  degrees  Kelvin  and  are 
totally  vaporized  (10:27).  The  X-ray  energy  is  completely 
absorbed  by  the  surrounding  air  at  distances  ranging  from 
ten  to  several  hundred  meters,  from  the  weapon,  creating  a 
spherical  fireball  that  initially  has  a  temperature  of 
several  thousand  degrees  Kelvin  (3). 

Because  the  fireball  is  hotter  and  less  dense  than  the 
surrounding  air,  a  large  pressure  differential  exists 
between  the  two  regions.  This  creates  a  strong  updraft 
leading  into  the  fireball.  If  the  burst  is  close  enough  to 
the  earth’s  surface,  large  amounts  of  dirt  and  dust  are 
sucked  by  the  updraft  up  into  the  fireball  and  are  vapor¬ 
ized  along  with  the  weapon  material.  The  fireball  also 
starts  to  rise  due  to  its  lower  density.  (10:35) 

t 

As  the  fireball  rises,  it  expands  and  cools  by  radiat¬ 
ing  energy  and  by  mixing  with  the  surrounding  air.  The 

dirt,  dust,  and  weapon  debris  condense  to  form  a  radioac- 

» 

tive  dust  cloud  that  has  now  become  toroidal  in  shape  due 
to  the  updrafts  (See  Fig.  1-1)  (10:28).  The  cloud  is  made 
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UPDRAFT  THROUGH 
CENTER  OF  TOROID 


Fig.  1-1.  Toroidal  Circulation  in  a 
Radioactive  Dust  Cloud  (10:29) 


up  of  particles  which  contain  radioactive  fission  products 
both  within  their  volumes  and  on  their  surfaces  (4).  Even¬ 
tually,  after  approximately  ten  minutes  or  so,  the  cloud 
cools  to  the  same  temperature  as  the  surrounding  air,  and 
is  said  to  be  "stabilized"  (10:32).  At  this  poin“,  the 
cloud  stops  rising,  although  it  will  continue  to  spread  in 
the  horizontal  direction  (10:32).  For  megaton  bursts,  the 
cloud  does  not  stabilize  until  it  has  reached  the  strato¬ 
sphere  . 

Under  the  influence  of  gravity,  the  radioactive  parti¬ 
cles  carried  aloft  and  formed  in  the  cloud  fall  to  the 
ground  at  rates  determined  by-  the  particle  size  (7,10). 
Where  the  particles  land  is  also  determined  by  particle 
size  and  by  the  local  winds.  The  deposition  of  these  par¬ 
ticles  on  the  ground  is  known  as  "fallout". 

The  amount  and  location  of  radioactive  fallout  are  of 
critical  importance  to  military  and  civilian  personnel 
involved  in  disaster  preparedness  or  other  similar  pro¬ 
grams.  Computer  codes  which  predict  fallout  have  existed 
since  the  late  1950's  (12:3).  Two  basic  categories  of 
fallout  models  exist  today:  numerical  disk  tossers  and 
smearing  models  (14:5.'. 

I 

The  numerical,  codes  model  the  cloud  formation  and 
particle  deposition  processes  in  much  greater  detail  than 
the  smearing  codes  do  and  are  able  to  handle  variable  wind 
fields  as  well  (22).  However,  they  are  complex  codes  that 
require  large  amounts  of  computer  time,  and  thus  are  not 
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suitable  for  studies  where  multiple  runs  of  the  code  are 
necessary  (23:1). 

The  smearing  codes  use  empirical  formulas  based  on 
fallout  data  from  nuclear  weapons  tests  (12:2).  They  are 
simple  and  fast,  but  cannot  handle  variable  winds,  assuming 
instead  a  constant  wind  in  one  direction  at  all  altitudes. 
This  is  a  critical  restriction  which  greatly  limits  the 
accuracy  of  these  models  (5:208),  Even  so,  smearing  codes 
are  the  codes  of  choice  for  large  studies  because  of  their 
fast  run  times.  Hopkins  developed  a  smearing  code  which 
treated  variable  winds  using  spectral  coefficients  (29). 
Hopkins'  code  proved  to  be  fgst  and  accurate  (15),  but  the 
code  was  developed  for  research,  is  not  documented,  and  is 
not  in  a  form  which  makes  it  easy  to  use. 

Problem  Statement 

The  goal  of  this  thesis  was  to  develop  and  document  a 
fast,  simple,  fallout  smear  code  that  uses  real  wind  fields 
to  transport  the  fallout  particles  to  the  ground  and  gives 
a  realistic  curved  hotline.  Hopkins'  method  of  spectral 
coefficients  is  used  to  model  the  wind  fields.  The  code  is 
intended  for  use  in  fallout  studies  at  the  Air  Force  Insti¬ 
tute  of  Technology,  Wright-Patterson  AFB,  Ohio. 

Scope 

The  computer  code  developed  in  this  thesis,  named 
« 

HYDRA  (for  Hotline  IM^ld-Dependent  Residual  Activity),  is 
limited  to  the  following  scenario: 


A 


1.  Only  one  weapon  can  be  detonated  at  any  given 
time. 

2.  The  weapon  is  detonated  on  the  earth's  sur¬ 
face.  Air  bursts  or  underground  bursts  can¬ 
not  be  modeled  using  this  code. 

3.  Only  local  fallout  is  modeled.  Local  fallout 
is  defined  as  fallout  which  reaches  the 
ground  in  24  hours  or  less. 

4.  Particle  transport  is  accomplished  by  gravity 
and  the  local  winds. 

5.  Dose  rate  contours  are  calculated  only  on  the 
earth's  surface.  No  attempt  is  made  to  calcu¬ 
late  the  activity  in  the  air  as  the  cloud 
descends  to  the  earth. 


Assumptions 


HYDRA  assumes  the  following: 

1.  The  fallout  particles  are  spherical  in  shape. 
Thus,  McDonald-Davies  fall  mechanics  can  be 
used  to  model  particle  fall. 

2.  Particles  do  not  begin  falling  out  of  the 
cloud  until  after  the  cloud  has  stabilized. 

3.  The  local  winds  change  gradually  over  dis¬ 
tance. 

4.  There  is  no  wind  component  in  the  vertical 
direction . 

5.  All  radioactivity  is  taken  into  the  cloud. 
There  is  no  stem  fallout. 

6.  There  is  no  washout  of  the  particles  from 
precipitation. 

General  Approach 

The  general  approach  taken  in  determining  the  fallout 
.  contours  consists  of  four  basic  steps.  First,  the  radioac- 

t 

tive  cloud  is  divided  into  a  number  of  particle  size  groups 
that  is  specified  by  the  code  user.  The  average  height  to 
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which  each  particle  size  group  is  injected  is  then  deter¬ 
mined  by  empirical  formulas  derived  by  Hopkins  from  numeric 
code  data  (14:14). 

Once  the  particle  injection  heights  have  been  calculat 
ed,  each  particle  size  group  is  transported  to  the  ground 
through  an  atmosphere  that  has  been  divided  into  ten  lay¬ 
ers,  with  the  center  of  the  top  layer  being  the  particle 
injection  height  for  that  particle  size.  Davies-McDonald 
fall  mechanics  (7,19)  are  used  to  determine  the  length  of 
time  it  takes  each  particle  group  to  reach  the  ground. 
Spectral  coefficients  are  used  to  calculate,  the  wind  fields 
in  each  atmospheric  layer  and  transport  the  particle  groups 
to  the  ground.  Connecting  the  ground  coordinates  for  each 
particle  group  yields  the  curved  hotline,  the  line  along 
which  the  dose  rate  will  be  a  maximum  for  a  given  distance 
from  the  burst  location. 

Once  the  hotline  has  been  located,  points  off  the  hot¬ 
line  are  found  for  user-specified  dose  rates  using  the  same 
empirical  formulas  found  in  the  smearing  codes. 

In  the  final  step,  the  hotline  points  and  dose  rate 
points  are  plotted  on  maps  to  give  dose  rate  contours  any¬ 
where  in  the  world.  DISSPLA  software  is  used  for  this  pur¬ 
pose.  The  progrart  which  takes  the  data  points  and  plots 
them  is  named  HYDMAP. 

User's  Guides  ’ 

The  development  of  user's  guides  for  HYDRA  and  its 
associated  codes  was  an  integral  part  of  this  project. 


These  guides  were  designed  to  be  used  with  the  CDC  Cyber 
750  computer  using  the  Network  Operating  System  (version 
2.4.2).  Three  user's  guides  have  been  written:  one  for 
HYDRA,  one  for  HYDMAP,  and  one  which  describes  how  to  cre¬ 
ate  a  file  containing  the  spectral  coefficients  from  raw 
wind  data. 

Sequence  of  Presentation 

The  remainder  of  this  report  is  divided  into  the 
following  sections: 

1.  Section  II  gives  brief  descriptions  of  some 
of  the  more  well-known  numeric  and  smearing 
codes . 

2.  Section  III  de'tails  the  program  code  theory 
and  development. 

3.  Section  IV  presents  the  program  evaluation 
and  discussion. 

4.  Section  V  presents  the  conclusions. 

5.  The  user's  guides  are  presented  in  Appendices 
D ,  E ,  and  F . 


II.  Current  Fallout  Prediction  Models 


Numeric  Codes 

Numeric  codes  employ  numeric  integration  to  determine 
ground  fallout.  The  radioactive  cloud,  the  atmosphere,  and 
time  are  divided  into  cells,  with  fallout  particles  being 
transported  to  the  ground  in  discrete  steps  (5:205).  The 
most  popular  and  well-known  of  the  numeric  codes  is  the 
DEfenne  Land  Fallout  Interpretive  £ode,  or  DELFIC.  DELFIC 
is  intended  for  use  in  research  and  to  be  a  standard  of  com 
parison  for  production  (analytic)  codes  (22:1).  It  "pre¬ 
dicts  local  fallout  from  nuclear  weapons  explosions  in  the 
yield  range  from  0.001  to  100,000  kilotons  over  a  range  of 
heights  of  burst  from  shallow  subsurface  to  fallout-safe 
airbursts"  (22).  Local  fallout  is  defined  in  the  DELFIC 
manual  as  "the  intensely  radioactive  material  which  falls 
to  the  ground  within  several  to  several  hundred  miles  of 
ground  zero"  (22:7). 

In  DELFIC,  the  radioactive  cloud  is  modeled  as  an 
"entraining  bubble  of  hot  gas"  by  dynamic  equations  using 
atmospheric  temperature,  pressure,  humidity,  and  wind  data 
(22:7).  After  the  cloud  has  stabilized,  it  is  divided  into 
several  cylindrical  shaped  clouds,  with  each  cloud  consist- 

i 

ing  of  particles  of  a  single  size  unique  to  that  cloud 
(22:27).  Each  of  the  particle  clouds  is  then  further  divid 
ed  in  the  vertical  direction  into  disks  (22:28). 

Once  the  different  disks  hove  been  identified,  each 
disk  is  transported  to  the  ground  using  variable  or  con- 
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stant  wind  fields.  The  atmosphere  is  divided  into  a  dis¬ 
crete  number  of  layers,  with  each  layer  subdivided  into 
rectangular  cells.  The  disks  are  then  transported  through 
this  layered  atmosphere,  their  horizontal  location  at  any 
time  being  determined  by  the  wind  vectors  associated  with 
each  rectangular  cell.  (27:10) 

The  rate  at  which  the  particles  fall  to  the  ground  is 
determined  by  the  particle  size.  DELFIC  assumes  all  parti¬ 
cles  are  spherical  and  that  the  particles  obey  the  Davies 
fall  mechanics  equations  (27:5).  The  particle  size  distri¬ 
bution  defaults  to  a  lognormal  distribution,  but  other  par¬ 
ticle  size  distributions  may  -be  specified  (22:15). 

Compared  to  the  smearing  codes,  DELFIC  models  the  phys 
ical  processes  of  a  nuclear  explosion  much  more  realistical 
ly  and  therefore  gives  more  realistic  fallout  contours. 
However,  as  mentioned  previously,  DELFIC  is  extremely  com¬ 
plex  and  requires  long  run  times,  and  therefore  is  not  suit 
ed  for  studies  where  the  code  must  be  run  repeatedly.  Also, 
the  contour  patterns  produced  by  DELFIC  have  been  observed 
to  break  up  at  distances  far  from  ground  zero  (2). 

Smearing  Codes  (WSEG-10  and  AFIT) 

Two  smearing  codes  will  be  reviewed  here:  1)  the  Weap¬ 
ons  System  Evaluation  Group  code  (WSEG-10),  and  2)  the  Air 
Force  Institute  of  Technology  code  (AFIT).  Both  use  empiri 
cal  formulas  derived  from  nuclear  weapons  test  data  (12:3). 
However,  the  AFIT  code  relies  less  on  empirical  formulas 
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than  does  WSEG-10.  In  the  AFIT  code,  empirical  relations 
are  used  mainly  for  describing  the  initial  location  of  the 
stabilized  cloud,  while  key  computations  a4re  based  on  physi 
cal  relations.  These  codes  are  known  as  "smear"  codes 
because  the  effect  achieved  by  the  dose  rate  formula  use'd 
is  to  smear  the  radioactive  cloud  along  the  ground  as  it 
travels  downwind  (5:208).  Smear  codes  have  very  fast  run 
times  and  no  contour  pattern  break-up,  but  the  accuracy  may 
be  less  than  that  achieved  by  the  disk  tosser  codes. 

The  dose  rate  formula,  the  fundamental  equation  of  the 
smearing  codes,  is  given  as  (5:207): 


^(x.y)  =  (K)(Yk)(FF)  /  f (x  ,  y , t ' )g( t '  )dt ' 


(II-l) 


where  D1(x,y)  is  the  unit  time  reference  dose  rate  in 
roentgens/hr , 


K  is  the  source  normalization  constant  in 
roentgen-kin^  per  hour-kiloton , 

Yk  is  the  weapon  yield  in  kilotons, 

FF  is  the  fission  fraction  (normally  token  to  be 
50%)  , 


f(x,y,t')  is  the  normalized  horizontal  distribution 
of  the  cloud  activity  per  unit  area  at 
time  t*  (units  of  km-^),  and 


g(t')  is  the  fraction  of  the  total  activity  that 

arrives  on  the  ground  per  unit  time  (units  of 
inverse  hours )  . 


The  term  Dj^x.y)  is  an  artificial  quantity  that  is  the 
dose  rate  "at  location  x  and  y  at  one  hour  if  all  of  the 
activity  which  will  eventually  land  at  x,y  is  assumed  to  be 
down  at  one  hour"  (5:207).  The  function  f(x,y,t')  is 
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assumed  to  be  Gaussian  in  both  the  x  and  y  directions  and 
is  given  by  the  formula  (5:207): 


1 


f (x,y ,t*  )  = 


V2ff  °x(t') 
1 


EXP 


-0.5 


x  -  vxt' 
ax(t') 


V2*  0  y ( t '  ) 


EXP 


\2 


-0.5 


y  -  vyt' 

Oy(t') 


(  1 1  —  2  ) 


where  ox(t')  is  the  standard  deviation  in  kilometers  of  the 
cloud  activity  in  the  x  direction, 

vx  is  the  x  component  of  the  wind  velocity  in  km/hr, 

Oy(t')  is  the  standard  deviation  in  kilometers  of  the 
cloud  activity  in  the  y  direction,  and 

Vy  is  the  y  component  6f  the  wind  velocity  in  km/hr. 


The  source  normalization  constant  is  a  function  of 
"the  radioactivity  produced  per  kiloton  of  fission  yield, 
the  energies  of  the  emitted  gamma  radiation,  the  properties 
of  the  air,  and  the  self-absorption  of  the  ground"  (5:207- 
208).  Various  values  for  the  source  normalization  constant 
have  been  used  over  the  years,  ranging  from  2654  to  6990 
r-km  /hr-kt  (8:24).  Davis  has  recently  calculated  a  lower 
bound  of  5729  r-km2/hr-kt  for  Nevada  Test  Site  soil  (8:23). 


WSEG-10.  WSEG-10  was  first  published  in  1959  and  has 
been  the  most  popular  smearing  fallout  prediction  code  in 
use  since  that  time.  The  code  models  fallout  from  nuclear 
explosions  with  yields  ranging  from  0.001  to  100  megatons, 
and  is  limited  to  surface  bursts.  (12:3-5) 

WSEG-10  models  the  radioactive  cloud  as  a  Gaussian 
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cylinder.  The  cloud  is  stabilized  and  the  center  of  the 
cloud  reaches  its  maximum  height  in  15  minutes  or  less 
(12:5).  The  maximum  cloud  center  height  is  given  by  the 
formula  (12:6): 


Hc  =  44.0  +  6 . lln( Ym) 

-  .205  | ln( Ym)  +  2 . 42 |  (ln(Ym)  +  2.42)  (II-3) 

where  Ym  is  the  weapon  yield  in  megatons  and  Hc  is  the  maxi¬ 
mum  cloud  center  height  in  kilofeet. 

The  g(t')  function  in  WSEG-10  is  an  empirical  fit 
to  weapons  test  data,  and  is  given  by  the  formula  (12:10): 

(F)  EX-P(t'/Tc) 

g(t’)  - ( H-4 ) 

(Tc)  F ( 1  +  l/n0) 


where  F  is  approximately  equal  to  1.0, 
T„  is  a  time  constant,  and 


n. 


=  1.5  -  . 25(Hc/60) . 


(11-5) 


WSEG-10  assumes  a  constant  wind  field  such  that  the 
wind  direction  becomes  the  x  direction  and  the  y  component 
of  the  wind  velocity  is  zero.  The  £(x,y,t')  function  then 
becomes : 


f  (  x  ,  y  ,  L  '  )  = 


>/  2*  *x(t') 


EXF 


-0.5 


\j 2*  .  o 


EXP 


x  ~  vxt' 

•x(t') 

,2  1 


v2  -i 


-0.5 


»y(f) 


(11-6) 


If  one  assumes  that  the  activity  is  constant  over  a  few 
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cloud  standard  deviations,  then  o„  ,(t’)  =  o  (t  ),  where  t 

X;y  x,yN  a 

is  the  time  of  arrival  of  the  cloud  center  line  on  the 
ground  for  a  given  particle  (4).  The  f(x,y,t*)  function 
can  then  be  split  into  an  f(x,t’)  function  and  an  f(y,t  ) 
function,  where: 


f  ( y ,  ta  ) 


V  2’ 


EXP 


°y  <  ca  > 


'  /  y  V 

-0.51 - 


( I 1-7 ) 


Since  f(y,t  )  is  no  longer  dependent  on  t1',  it  can  be  moved 
outside  the  integral  in  Equation  I I — 1 ,  which  then  becomes: 


Djl  (  x  ,  y  )  =  (K)(Yk)(FF)£(y,t  ) 


/OO 

f  (x, 

n 


t ' )g( t  * )dt  *  ( II-8) 


If  it  is  assumed  that  the  g(t')  function  is  a  slowly 
changing  function,  then  it  can  be  represented  in  a  Taylor 

expansion  about  the  time  of  arrival  (4).  If  one  changes 
the  lower  limit  of  integration  in  Equation  11-8  from  0  to 
~  <o  ,  the  integral  may  be  solved  analytically  with  the  solu¬ 
tion  being  f  4 ) : 

J f(x,t')g(t')dt'  *  g(ta)/vx 

-CO 

Equation  11-8  then  becomes  an  analytical  expression  easily 
solved  in  a  computer  code  (14:8): 


(K)(Yk)(FF)g(ta) 

D^x.y)  =  - - -  EXP 


\j  2* 


ffy(ta)vx 


2  I 


-0.5 


a  (  t  ) 
y  x  a  7 


(11-9) 


The  activity-size  distribution  for  the  WSEG-10  model 
which  led  to  the  empirical  expression  for  g(t)  was  (12:7): 


13 


1 


EXP 


( 11-10) 


A(r) 


yj 2*  0  r 


1  n C rm  )  -  ln(r)vi 


0 


where 


r 

r 

0 


is  the  median  radius  (=  44  microns), 
is  the  particle  radius  in  microns,  and 
=  0.690. 


While  fast  and  simple,  WSEG-10  has  many  limitations 
which  adversely  affect  its  accuracy.  These  limitations 
were  summarized  by  Bridgman  and  Bigelow  (5): 

1.  Restriction  to  a  constant  wind  field  gives  poor 
results  if  winds  are  varying. 

2.  The  form  of  the  g(t')  function  predicts  activi¬ 
ty  on  the  ground  at  zero  time,  which  is  physi¬ 
cally  unrealistic. 

3.  The  activity-size  distribution  assumes  all 
activity  is  evenly  distributed  throughout  the 
particle  volume.  In  fact,  it  is  fractionated, 
with  part  of  the  activity  distributed  in  the 
volume  and  the  rest  deposited  on  the  surface. 

4.  The  activity-size  distribution  predicts  too 
high  an  activity  level  downwind  by  putting  too 
much  of  the  activity  in  particles  with  radii  of 
10  microns  or  more  (the  median  radius  is  too 
large) . 

5.  For  low  yield  weapons,  the  downwind  activity 
close  to  ground  zero  is  underestimated  because 
the  activity-size  distribution  is  too  narrow 

(  0  is  too  small)  . 

6.  The  particle  size  cannot  be  varied  (2).  This 

is  especially  limiting  since  the  g(t)  function 
(as  will  be  shown  later)  is  dependent  on  the 
particle  size.  , 


A  FIT .  The  AFIT  code,  developed  at  the  Air  Force  Insti¬ 
tute  of  Technology  also  uses  Equation  II-9,  but  contains 
significant  improvements  which  answer  criticisms  2-6  listed 
in  the  previous  section.  The  improvements  were  made  in  the 
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activity-size  distribution  and  in  the  g(t')  function. 

The  AFIT  normalized  activity-size  distribution  is 
actually  the  combination  of  two  lognormal  functions  and  is 
given  by  the  formula  (5:211): 


A  ( r  ) 


f 


v 


\] 2*  0  r 


+ 


(1-f  ) 
v  v 

yjl*  0  r 


(11-11) 


where  fv  is  the  fraction  of  the  activity  that  is  dis¬ 
tributed  throughout  the  particle  volume, 

a  2  =  ^ 11  (  r  m)  +  2  0  j  , 

a3=  ln(rm)  +  3  P  ,  ; 

r  .is  the  median  size  of  the  number-size  distribu- 
m  .  . 
t  x  o  n ,  and 

0  is  the  logarithmic  slope  of  the  number-size  dis¬ 
tribution. 


In  the  AFIT  A(r)  distribution,  rm  and  0  can  be  varied, 
whereas  in  the  WSEG-10  code,  these  values  are  fixed. 

The  AFIT  g(t)  function  (dropping  the  prime)  is  based 
on  McDonald-Da vies  fall  mechanics  (7,10),  and  is  not  an 
empirical  fit  to  weapons  test  data  as  is  the  WSEG-10  g(t) 
function.  Bridgman  reasoned  that  g(t)  and  A(r)  are  both 
particle  distribution  functions  and  thus  must  be  related 
according  to  the  equation  (5: 210) :* 

g(t)dt  =  A(r)dr 

or  ,g(t)  =  A'r)dr/dt  (11-12) 


where  dr/dt  is  the  rate  of  ch  nge  in  the  size  of  the  par¬ 
ticles  hitting  the  ground  with  respect  to  the  time  of  arri 


val.  If  A(r)  is  known,  all  one  must  do  is  find  dr/dt  to 
determine  g(t).  In  one  version  of  the  AFIT  code,  the  quan¬ 
tity  dr/dt  is  determined  using  sets  of  coefficients  pub¬ 
lished  by  Colarco  (6).  These  coefficients  resulted  from  an 
empirical  fit  to  results  calculated  from  the  McDonald- 
Davies  fall  mechanics  equations  described  in  the  next  chap¬ 
ter.  The  coefficients  are  dependent  on  the  particle  injec¬ 
tion  height  (6:7).  The  AFIT  g(t)  function  has  a  value  of 
zero  when  the  arrival  time  is  zero,  thus  answering  criti¬ 
cism  #2  of  the  previous  section. 
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Overview 

The  program  developed  for  this  thesis,  named  HYDRA 
(Hotline  Yield-Dependent  Residual  Activity),  is  essentially 
identical  to  Hopkin's  variable  wind  code  (14).  The  fallout 
cloud  is  divided  into  a  user-specf ied  number  of  particle 
size  groups  which  are  gravity-sorted  within  the  cloud,  the 
heavier  particles  residing  in  the  lower  part  of  the  cloud. 
Each  particle  size  group  is  then  transported  separately  to 
the  ground  through  a  discretely  layered  atmosphere  using 
McDonald-Davies  fall  mechanics  (7,19)  and  real  winds.  Wind 
field  vectors  are  determined  using  the  method  of  spectral 
coefficients  at  specified  points  along  a  given  particle 
trajectory  (14,15). 

The  connected  ground  coordinates  for  each  particle 
size  group  form  the  fallout  hotline.  Along  this  hotline, 
the  dose  rate  will  be  a  maximum  for  any  specified  distance 
from  the  weapon  detonation  point.  The  dose  rate  for  any 
point  on  or  off  the  hotline  is  calculated  using  Equation  II- 
1.  Dose  rate  contours  are  then  overlayed  on  maps  to  indi¬ 
cate  local  fallout  for  any  point  on  the  globe  excluding  the 
poles. 

The  program  code  theory  and  program  development  is 

divided  into  the  following  sections: 

« 

1.  Particle  size  selection. 

2.  Cloud  rise  modeling. 
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3.  Particle  fall  mechanics. 

4.  Activity-size  distribution. 

5.  Wind  vector  calculations. 

6.  Wind  shear  calculations. 

7.  Calculation  of  g(t). 

8.  Dose  rate  calculations. 

9.  Dose  rate  contour  mapping. 

Each  topic  is  discussed  in  one  of  the  following  sections. 

Particle  Size  Selection 

HYDRA  is  intended  for  local  fallout  studies,  where 
local  fallout  is  the  fallout  that  reaches  the  ground  in  24 
hours  or  less.  Having  the  user  choose  the  particle  sizes 
which  all  fall  within  24  hours  may  lead  to  difficulties, 
since  the  arrival  time  for  a  given  particle  size  depends  on 
the  weapon  yield.  Therefore,  a  routine  was  developed  that 
automatically  selects  the  particle  sizes  once  the  number  oi 
particles  ha3  been  specified  by  the  user.  This  routine 
first  finds  the  radius  of  the  particle  that  will  lsnd  in 
approximately  24  hours  (hereafter  referred  to  as  the  criti¬ 
cal  radius),  and  calculates  the  corresponding  value  of  the 
cumulative  lognormal  activity-size  distribution  (See  Fig. 
III-l).  The  range  between  this  corresponding  value  and  the 
highest  possible  value  of  1.0  is  then  divided  into  a  number 
of  areas  equal  to  the  number  of  particles,  with  each  area 
being  represented  by  its  corresponding  particle  radius. 

The  critical  radius  is  found  using  McDonald-Davies 
fall  mechanics.  The  particle  injection  height  is  deter- 


I’iy.  111-1.  Cumulative  Activity-Si2e  Distribution 


mined  using  Equation  II-3  from  WSEG-10.  The  particle  falls 
from  this  height  to  the  ground  through  a  single  layered 
atmosphere.  The  average  fall  velocity  for  the  particle  is 
assumed  to  be  the  velocity  calculated  at  the  altitude  mid¬ 
way  between  the  injection  height  and  sea  level.  The  time 
of  arrival  to  the  ground  is  then  simply  the  injection 
height  divided  by  the  fall  velocity.  An  initial  particle 
radius  of  5.0  microns  is  chosen,  and  its  time  of  arrival  is 
calculated.  The  radius  of  the  particle  is  increased  in 
increments  of  1.0  microns  in  a  simple  iterative  process 
untiL  the  arrival  time  is  less  than  24  hours.  One  micron 
is  then  added  to  the  particle  radius,  and  this  new  particle 
radius  is  the  critical  radius. 

Once  the  critical  radius  has  been  determined,  the  cor¬ 
responding  value  of  the  cumulative  activity-size  distribu¬ 
tion  is  calculated.  The  activity-size  distribution  used 
here  is  Freiling's  approximation  of  Equation  X I — 1 1  (4): 


(1II-1) 


x  =  (ln(r)  -  a  2.5)  /  0 


( I I I— 3 ) 
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then  Equation  1 1 1  —  2  becomes: 

Cum.  A(x)=/A(x')dx'  ( I I I — 4 ) 

J  —OO 

where  A(x)  is  a  normal  function.  Equation  III-4  is  solved 
using  an  IMSL  Fortran  subroutine  to  find  the  cumulative 
distribution  value  corresponding  to  the  critical  particle. 

Once  the  cumulative  distribution  v.alue  of  the  critical 
radius  has  been  determined,  the  range  between  this  distribu¬ 
tion  value  and  the  highest  possible  distribution  value 
(1.0)  is  divided  into  equal  segments  whose  number  is  equiva¬ 
lent  to  the  number  of  particles.  Then,  the  cumulative  func¬ 
tion  values  representing  each  segment  are  fed  into  another 
IMSL  subroutine  to  find  the  corresponding  particle  radii. 

Cloud  Rise  Modeling 

The  cloud  rise  model  used  in  HYDRA  is  similar  to  the 
one  found  in  WSEG-10  and  the  AFIT  model.  The  cloud  is 
assumed  to  stabilize  quickly  and  particle  fallout  does  not 
begin  until  after  the  cloud  has  stabilized.  However,  where¬ 
as  WSEG-10  and  earlier  versions  of  the  AFIT  code  use  a  sin¬ 
gle  particle  injection  height  for  all  particle  sizes,  HYDRA 
uses  a  different  injection  height  for  each  particle  size 
group.  Each  height  represents  the(  center  of  a  wafer  cloud 
consisting  of  monosized  particles  (14:11). 

Hopkins  discovered  that  using  a  single  injection 
height  for  all  particle  sizes  "severely  restricts  hotline 
curvature"  (14:31),  giving  a  false  hotline.  His  solution 


to  the  problem  was  to  derive  an  empirical  formula  for  par¬ 
ticle  injection  heights  based  on  a  least  squares  fit  to 
DELFIC  data  for  yields  ranging  from  1.0  kilotons  to  15.0 
megatons  (14:11).  The  formula  is  (14:14): 

H  =  (SLOPE) ( D  )  +  INTERCEPT  (III-5) 

c  p 

where  Ilc  is  the  average  height  in  meters  of  the  wafer  cen¬ 
ter  in  the  stabilized  cloud, 

Dp  is  the  particle  diameter  in  microns, 

SLOPE  =  -EXP( 1.574  -  .011971nY  +  .03636(lnY)2 

-  . 004 1 ( In Y )3  +  .0001965(lnY)4) 
in  meters/ micrometer, 

INTERCEPT  =  EXP( 7 . 889  +  .341nY  +  . 001 226 ( InY )2 

-  . 005227 (InY)3  +  . 0004 1 7 ( InY )4 ) 

is  the  altitude  in  meters  for  a  particle  of 
zero  diameter  on  a  linear  fit,  and 

Y  is  the  weapon  yield  in  kilotons. 

The  particle  injection  height  is  determined  for  each 

particle  in  subroutine  ARRTIM. 

Particle  Fall  Mechanics 

The  particle  fall  mechanics  are  modeled  using  the 
McDonald -Da vies  equations  and  a  discretely  layered  atmos¬ 
phere.  All  particles  are  assumed  to  be  spherical.  First, 
the  atmosphere  is' divided  into  ten  layers  as  in  llopkin's 
variable  wind  model.  Hopkins  found  that  while  the  optimum 
number  of  layers’was  24,  using  10  layers  resulted  in  only  a 
5.0  percent  decrease  in  accuracy  while  reducing  the  comput¬ 
er  calculation  time  by  a  factor  of  two  (24:29).  The  center 
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of  the  highest  layer  is  the  particle  injection  height,  so 
each  particle  size  has  its  own  set  of  ten  layers.  This  is 
different  from  Hopkins'  model,  where  one  set  of  layers  is 
defined  for  all  particle  sizes,  with  the  center  of  the  high¬ 
est  layer  being  the  injection  height  of  the  smallest  parti¬ 
cle  . 

For  the  top  and  bottom  of  each  layer,  the  air  density, 
dynamic  viscosity,  and  particle  fall  velocity  are  calculat¬ 
ed.  Air  density  and  dynamic  viscosity  are  determined  using 
empirical  formulas  of  the  U.S.  Standard  Atmosphere  (21). 

The  McDonald-Davies  equations  are  used  to  determine  the 
particle  fall  velocity.  The  force  balance  equation  for 
falling  particles  in  air  is  (19:463): 

.5  PgV^C^jrr2  =  (4/3)7rr3  pfg  (III-6) 

where  pa  is  the  air  density,  is  the  drag  coefficient,  vz 
is  the  fall  velocity,  r  is  the  particle  radius,  and  the 

O 

quantity  .5pav§  is  the  dynamic  pressure.  For  a  given  par¬ 
ticle  size  and  altitude,  there  are  two  unknowns  in  this 
equation:  vz  and  Cd  .  The  Reynolds  number  for  falling 
spheres  is  given  by  (19:463): 

R  =  2y  p  r/  '  (III-7) 

L  a  x 

where  q  is  the  dynamic  viscosity.  Substituting  Equation 
1 1 1  —  7  into  Equation  1 1 1 —6  yields  the  following  relationship 

t 

(5:212) : 


R2Cd  =  32papfgr3/(3>,2) 


(III-8) 


Thus,  the  quantity  R2Cd  can  be  easily  determined  for  a  giv¬ 
en  particle  size  and  altitude. 

Once  the  quantity  R2 Cd  is  known,  the  Reynolds  number  R 
can  be  determined  using  empirical  relations  developed  by 
Davies  (5:212): 

R  =  R2Cd/24  -  2.3363E-04(R2Cd)2  +  2 . 0154E-06 (R2 Cd  )  3 
-  6 . 9105E-09 (R2  Cd ) 4 
for  R2Cd  <  100.0,  or 

log(R)  =  -1.29536  +  .9861og(R2Cd)  -  .  046677  [  log  (R2  Cd)  ]  2 
+  .0011235[log(R2Cd )]3 

for  R2Cd  >  100.0  ( III-9) 

Equation  III-7  can  then  be  used  to  find  the  fall  velocity 
at  a  given  layer  boundary.  This  value  is  then  corrected 
for  drag  slip  at  high  altitudes  by  multiplying  it  by  a  cor¬ 
rection  factor  taken  from  DELFIC  (5:212): 

CF  =  1  +  1 . 165E-07/ r  p  a  (III-10) 

where  r  is  the  particle  radius  in  meters  and  pa  is  the  air 
density  in  kg/m3. 

The  particle  fall  velocity  through  a  given  layer  is 

assumed  to  be  the  average  of  the  velocities  at  the  layer 

boundaries.  The  time  it  takes  the  particle  to  fall  through 

the  layer  is  then  simply  the  layer  thickness  divided  by  the 

average  fall  velocity.  Summing  the  individual  fall  times 

* 

for  each  layer  gives  the  time  of  arrival  on  the  ground  for 
a  particle. 
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Air  density  and  dynamic  viscosity  at  the  layer  boundar¬ 
ies  are  calculated  in  subroutine  STDATM.  The  fall  velocity 
at  the  layer  boundaries  is  calculated  in  subroutine  VFALL. 
Both  subroutines  are  called  by  subroutine  ARRTIM.  Average 
fall  velocities,  fall  times,  and  the  time  of  arrival  are 
calculated  in  subroutine  ARRTIM. 

Activity-Size  Distribution 

The  activity-size  distribution  used  in  HYDRA  is  the 

same  one  used  in  the  AFIT  code,  and  is  given  by  Equation  II- 

11.  It  is  calculated  in  subroutine  ACTVTY.  The  parameters 

r  and  0  are  input  variables, 
m 

The  g(t)  Function 

The  g(t)  function  in  HYDRA,  as  in  the  AFIT  code,  is 
given  by  the  following  equation  (5:210): 

g(t)  =  A(r)dr/dt  (III-U) 

However,  a  different  method  is  used  to  calculate  dr/dt. 

Where  the  AFIT  code  employed  the  Colarco  coefficients 
to  calculate  dr/dt,  HYDRA  uses  a  finite  difference  method. 
First,  the  given  particle  size  is  multiplied  by  0.99  and 
1.01  to  yield  two  particles  sizes  which  bound  the  given  par¬ 
ticle  size.  Then,  the  arrival  times  for  these  "bounding" 
particles  are  calculated  using  the  methods  described  in  the 

above  sections.  The  value  for  dr/dt  is  easily  calculated 

* 

using  the  formula: 

dr/dt  -  ABS  ( R2  -  R1)/(T2  -  Tl)  (111-12) 
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where  R2  =  1.01  X  (particle  radius), 

R1  =  0.99  X  (particle  radius), 

T2  is  the  arrival  time  for  the  particle  of  radius  R2, 
and  T1  is  the  arrival  time  for  the  particle  of  radius  Rl. 

Colarco’s  coefficients  were  not  used  because  Bridgman 
reported  deviations  of  10  to  13  percent  in  the  arrival 
times  for  certain  radii  when  compared  to  pure  fall  mechan¬ 
ics  theory.  (2) 

The  quantity  dr/dt  is  calculated  in  the  main  program 
of  HYDRA.  The  g(t)  function  is  calculated  in  subroutine 
ACTVTY . 

Wind  Shear  Calculations 

Wind  shear  calculations  in  HYDRA  are  identical  to 
those  in  Hopkins'  model.  Wind  shear  is  the  "local  varia¬ 
tion  of  the  wind  vector  or  any  of  its  components  in  a  given 
direction"  (17:631).  For  a  discretely  layered  atmosphere, 
the  shear  of  a  wind  component  in  any  given  layer  is  the 
difference  in  the  wind  components  at  the  top  and  bottom  of 
the  layer  divided  by  the  layer  thickness.  However,  in 
HYDRA  the  shear  in  a  given  layer  is  weighted  by  the  resi¬ 
dence  time  of  the  particle  in  the  layer  such  that  the  total 
shear  experienced  by  a  particle  can  be  represented  by  the 
formula  (15:36):  . 

SHRTOT  =’[E(SHR2(i)  x  TFALL(i)/TARRIV)J  *  (III-13) 

where  SHRTOT  is  the  total  root-mean-square  shear  in  the  X 
or  Y  direction, 


SHR(i)  is  the  shear  experienced  by  a  particle  in  the 
ith  layer  in  the  X  or  Y  direction, 

TFALL(i)  is  the  residence  time  of  the  particle  in  the 
ith  layer, 

TARRIV  is  the  arrival  time  of  the  particle  on  the 
ground , 

and  TFALL(i)/TARRIV  is  the  weighting  factor  for  the  ith 
layer . 

Wind  shear  will  spread  the  fallout  contour  out  lateral 
ly  from  the  hotline  (18:1).  All  wind  shear  calculations 
are  performed  in  the  main  program  of  HYDRA. 

Variable  Wind  Calculations 

The  distance  in  the  X  and  Y  directions  that  a  particle 
travels  as  it  falls  through  a  given  layer  is  assumed  to  be 
the  particle  fall  time  for  that  layer  multiplied  by  the 
average  wind  velocity  components  encountered  in  the  layer. 
The  distances  computed  for  each  layer  are  then  summed  to 
find  the  total  horizontal  transport  experienced  by  the 
particle  as  it  falls  to  the  ground. 

The  wind  velocity  components  may  be  determined  in  one 
of  two  ways.  The  first  method,  employed  in  DELFIC,  uses  a 
three  dimensional  grid  and  assigns  wind  component  values  to 
the  grid  points.  Values  for  locations  not  on  the  grid 
points  are  linearly  interpolated.  'The  second  method,  used 

i 

in  Hopkins'  code  pnd  in  HYDRA,  represents  the  wind  compo¬ 
nents  as  the  truncated  sura  of  orthogonal  functions  in  spher 

* 

ical  coordinates,  and  is  known  as  a  spectral  representa¬ 
tion.  To  use  this  method,  the  complex  coefficients  of  the 
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functions  must  be  determined.  (11:214) 

The  main  reason  that  spectral  representation  was  chos¬ 
en  over  the  grid  method  is  that  an  evenly  spaced  gridded 
wind  mesh  fine  enough  for  linear  interpolation  is  not 
available  on  a  global  scale.  Spectral  representation  has 
an  additional  advantage  in  that  there  is  no  need  to  map  the 
globe  (1:399). 

The  complex  coefficients  are  derived  from  wind  data 
collected  at  1415  stations  around  the  northern  hemisphere. 
Balloons  with  sensors  are  released  at  these  stations  twice 
a  day  at  midnight  and  noon  Greenwich  Mean  Time  (GMT).  The 
lateral  displacement  of  each  .balloon  is  measured  at  certain 
altitudes,  and  from  these  lateral  displacements  the  wind 
data  is  calculated.  The  wind  data  is  fed  into  the  program 
AFGL  which  outputs  the  complex  wind  coefficients.  AFGL  is 
composed  of  subroutines  taken  from  a  weather  forcasting 
model  used  at  the  Air  Force  Geophysics  Laboratory,  Ilanscom 
AFB.  The  complex  wind  coefficients  are  contained  in  the 
file  SPECOEF,  which  is  an  input  file  for  HYDRA.  These  wind 
coefficients  are  divided  into  12  sets  or  "spectral  layers", 
one  set  for  each  of  the  12  altitudes  at  which  the  wind  vel¬ 
ocity  is  determined  (See  Fig.  II1-2).  (24:2,3) 

The  vorticity,  divergence,  continuity,  adiabatic  ther¬ 
modynamic,  and  hydrostatic  equations  are  the  set  of  equa¬ 
tions  that  describe  the  adiabatic,  frictionless  flow  of  the 
atmosphere  (20:737).  The  vorticity  equation  is  the  curl  of 
the  vector  equation  of  motion,  while  the  divergence  eqea- 
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tion  results  from  taking  the  divergence  of  the  equation  of 
motion  (17:174,615). 

The  above  equations  all  contain  the  velocity  in  vector 
form  as  a  variable.  The  Helmholtz  theorem  states  that  the 
horizontal  wind  vector  can  be  expressed  as  the  sum  of  irro- 
tational  and  non-divergent  fields  according  to  the  formula 
(20:737): 

V  =  +  Vx  (III— 14) 

with 

Y.  =  1c  X  ;  Vx=  VX  (HI-15) 

where  i  is  the  stream  function, 

X  is  the  velocity  potential, 

1c  is  the  vertical  unit  operator,  and 
V  is  the  horizontal  gradient  operator. 

The  stream  function  is  a  scalar  function,  and  for  non- 
divergent  flow  is  defined  such  that  (17:548): 

u  =  —  l  ;  v  =  (III— 16) 

where  u  is  the  wind  component  in  the  X  or  longitudinal 
direction,  and  v  is  the  wind  component  in  the  Y  or  latitudi¬ 
nal  direction.  The  velocity  potential  is  also  a  scalar 
function,  and  for  irrotational  flow  is  defined  such  that 
(17:608): 

u  -  -  dx  /dx  ;  v  »  —  ax  /8y  (III-17) 


30 


The  quantities  u  and  v  are  the  ones  needed  to  compute 
the  horizontal  particle  transport  in  HYDRA.  However,  the 
attempt  to  represent  these  quantities  as  spectral  sums  is 
an  awkward  exercise  which  can  be  avoided  if  the  following 
substitutions  are  made  (24:8): 


U  =  u(sin  0  )  ;  V  =  v(sin9)  (III  — 18) 

»  l  , 

where  0  is  the  colatitude. 

Now,  the  task  is  to  represent  the  quantities  U  and  V 
as  sums  of  the  orthogonal  functions  Y^ ,  where  Y^  are  func¬ 
tions  of  the  longitude  X  and  colatitude  0.  These  func¬ 
tions  Y^  are  solutions  of  the  equation  (11:215): 

a2  V  V  +  n(n+l)Y™  =  0  (111-19) 

where  a  is  the  radius  of  the  earth.  The  Y™  functions  are 
of  the  form  (11:215): 


yJJ  =  EXP( im X  )  X  pjj  (111-20) 


where  1^  are  polynomials  of  order  m  and  degree  n.  Substi¬ 
tuting  Equation  III-20  into  Equation  1 1 1  —  1 9  yields 
(11:215) : 


d2Pj  dPn' 

- - —  +  cot©  - 


d0 


•  d  t 


m 


n(n+l)  - 


sin20 


P*  =  0 


(111-21) 


where  t  represents  time. 

Equation  J 11—21  is  a  variation  of  Legendre's  equation, 
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and  its  solutions  are  associated  Legendre  polynomials 

(11:215).  The  solutions  Pm  ,  when  normalized,  are  defined 

n 


as  (29:24): 


pm 

n 


(-D 


m+n 


2nn  ! 


( 2n+l ) ( n-m) ! 


2 ( n+m) ! 


(1-X2)”/2 


X 


dm+n 

dxm+n 


(l-x2)n 


(III-22) 


where  x  =  cos  0  .  In  HYDRA,  the  functions  Pm  are  calculated 

n 

using  the  recursion  relation  (29:24): 


(m  p  m 
n+1  n+1 


Cx)  +  ( 


m  pm 
n  n-1 


(x) 


where  (29:24): 


m 

< 

n 


4n^  -  1 


( III-23) 


( III-24) 


The  values  for  <mare  evaluated  in  subroutine  EPSLON,  while 

n 

the  Pm  functions  are  evaluated  in  subroutine  PLN3.  These 
n 

subroutines  were  taken  directly  from  Ilopkin's  code. 

A  real,  smooth  function  F(  X  ,0)  that  is  defined  for 
all  values  of  X  and  0  over  a  sphere  can  be  expressed  as 
the  sum  of  a  series  of  orthonormal,  functions  as  shown 
(13:146,9:263):  ' 

CO  00 

F(  X  ,0)  =  £  T  FmYm  ( II 1-25) 

*  m=--oo  n=jm|  n  n 

where  FJJ  are  complex  expansion  coefficients  and  Y™  is 
defined  by  Equation  III-20.  However,  Horn  the  definition 


'  t  *_■  .m_?  ry  *,»  t  ,*_■  /■_»  .. ■*  „» .  * 
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of  P™  in  Equation  III-22,  is  zero  for  |m|  >  n  (20:737). 
Thus,  the  summation  in  Equation  III-25  can  be  limited  as 
follows : 


n  oo 

F(  X  ,0)  =  £  £  ,  F"^  ( 1 1 1  —  2 6 ) 

m=-n  n^Jm|  n  n 

Since  the  summation  limits  must  be  finite  for  the  program 
to  compute  values  for  the  variables,  the  limits  are  set  to 
a  number  J  such  that  (29:25): 


F(  A  ,8) 


pm  y  ^ 

n  n 


(III  —  27) 


The  truncated  summation  in  Equation  III-27  is  known  as  a 
rhomboidal  truncation  because  plotting  the  numbers  n  versus 
m  results  in  a  rhomboid  shape  (See  Fig.  III-3).  Rhomboidal 
truncation  is  used  in  this  situation  because  it  allows  the 
wind  vectors  to  be  resolved  to  the  same  degree  in  both  the  X 
and  0  directions  (29:2).  In  HYDRA,  J  is  set  to  30,  the 
same  truncation  limit  used  in  atmospheric  models  at  the 
National  Meteorological  Center. 

Using  Equation  111-27,  U  and  V  can  be  represented  by 
the  following  equation  (28:1284): 

J  | mj+J+1 

(U,V)  =  £  .  £  (Um,Vm)Ym  (111-28) 

m=-J  n=H  n  n  n 

where  the  summation  has  been  modified  slightly  to  keep  the 
U  and  V  terms  compatible  with  other  a tjnospher ic  variables 
using  the  truncation  limits  of  Equation  1 1 1  —  2  7  (28:1284). 

t 

One  further  simplication  of  the  summation  limits  may 
be  made.  The  summation  over  m  may  be  broken  down  as  fol- 
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Fig.  III-3.  Spectral  Domain 
Using  a  Rhomboidal  Truncation  (24:13) 
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lows  (15:122): 

(ni-29) 

m=-J  m=-J  m=0  m=l 

The  summation  over  the  negative  values  of  m  can  be  rewrit¬ 
ten  as  (15:122): 

-1  | m| +J+1  J  m+J+1 

E  E  <U"'V")Yn  =  E  E  dC'C^n”  (HI-30) 

m=-J  n=|m|  m=l  n=m  '  ' • 

However,  since  U  and  V  are  real,  physical  quantities,  the 

spectral  coefficients  (U  m,V  m )  can  be  represented  as 

n  n 

(29:25) : 


(U~m,V~m)  =  (-l)m  (um,vm)* 


( 1 1 1-31 ) 


where  (Um,Vm)*  are  the  complex  conjugates  of  the  spectral 
coefficients.  Likewise  (15:123): 


y-m  _  /  | s-m (Ym)+ 

n  v  '  s  w' 


(III-32) 


where  (Y™)*  is  the  complex  conjugate  of  Y™.  The  summation 
over  the  negative  values  is  then  (15:123): 

J  m+J+1 

E  E  <Unm  -V;m)Y;m  = 

m=l  n=m 

J  m+J +1 

E  E  (CO*  (Y®)*  (111-33) 

m=l  n=m  u 


4’V 

L  « • 

vjy 


Since  U  and  V  are  real  quantities ,,  only  the  real  part  of 
Equation  111-33  is  used  (15:123): 

J  m+J+1 

Re  E  E  (Um,Vm)*  (Y®)*  = 
m=l  •  n=m  n  n  n 


J 

Re  E 


m=l  n=m 


m+J+1 

E  (Um,Vm)Ym 

*-•  v  n  ’  n  '  n 


(111-34) 
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This  is  equivalent  to  the  real  part  of  Equation  1 1 1  —  2 8  when 
summed  over  positive  values  of  m.  Therefore,  the  final 

form  for  U  and  V  can  be  given  by  the  formula  (15:124): 

J  m+ J +1 

01. V)  =  Re  Y,  H  A(U£,V£)Y£  (111-35) 

m=0  n=m 

where  A  =  1  if  m  =  0  and  A  =  2  for  nonzero  values  of  m. 
Equation  III-35  is  the  formula  used  in  HYDRA  to  calculate 
the  wind  components  at  a  given  spectral  layer  altitude,  and 
is  found  in  subroutine  UVCOMP,  which  was  taken  directly 
from  Hopkin's  code.  Wind  velocity  components  at  altitudes 
other  than  the  spectral  layer  altitudes  are  linear  interpo¬ 
lations  of  the  wind  components  of  the  spectral  layers  which 
bound  that  altitude. 


Dose  Rate  Calculations 


Dose  rate  calculations  are  performed  assuming  the 
cloud  distribution  is  Gaussian  in  the  horizontal  and  verti¬ 
cal  directions.  HYDRA  uses  Equation  I I — 1  to  determine  the 
dose  rate  both  on  and  off  the  hotline,  but  the  final  form 
is  different  from  that  used  in  WSEG-10  and  the  AFIT  model 
because  of  the  curved  hotline.  The  final  form  is  (14:82): 

2.,  2 


D1(x,y)  = 


(K)(Yk)(FF)g(g 


yj  2ti 


2  v  2 
*yX 


*xY 


EXP 


-.5 


(xY  -  yX) 


2y2  . 
cr  yX  + 


2„  2 


(111-36) 


where  (X,Y)  is  the  particle  arrival  point  on  the  hotline 
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and  (x,y)  is  a  point  off  the  hotline.  This  form  arises 
mainly  because  the  vx  t '  and  Vy  t '  terms,  where  vx  and  Vy  are 
constants,  are  replaced  with  the  integrals  (14:73): 

V  =>  /vxdt' 

°rv 

V  =>  ivydt' 

For  a  detailed  derivation  of  Equation  III.-36,  see  Reference 
14,  p73-82.  A  comparison  of  Equation  III-36  with  Equation 
II-9  will  show  that  the  term  raised  to  the  -  power  in  Equa 
tion  III-36  is  functionally  equivalent  to  the  ffy(t')vx  term 
in  Equation  II-9. 

Equation  III-36  models  t-he  grounded  activity  as  a  Gaus 
sian  distribution  centered  about  the  hotline.  St  Ledger 
showed  that  this  Gaussian  distribution  lies  on  a  straight 
line  that  is  perpendicular  to  a  line  extending  from  ground 
zero  to  a  particle  arrival  point  on  the  hotline  (30:32- 
35).  This  straight  line  passes  through  the  arrival  point, 
and  its  equation  is  given  by: 

y  =  ( -X/Y ) x  +  Y  +  (X2/Y)  (III-37) 

In  calculating  contour  points  off  the  hotline,  the  dose 
rate  D^(x,y)  is  known  and  it  is  the  co.ordinates  x  and  y 

i 

which  must  be  found.  Substituting  Equation  III-37  into 

Equation  1 1 1  —  3  6  yields  the  coordinate  as  a  function  of 

L) .  ,  X,  Y,  Y.  ,  K,  a  ,  a  ,  and  t  .  The  coordinate  y  can 
i  k  i  x  y  & 

then  be  found  using  Equation  III-37. 

The  terms  a  (t  )  and  a  ( t  )  ,  the  cloud  deviations  in 


the  x  and  y  directions,  arise  from  the  toroidal  motion  in 
the  cloud  and  from  wind  shear  (4).  The  equation  for  these 
terms  is  ( 4 ) : 


a 


2 


(ta)’- 


+ 


(III-38) 


where  S  is  the  shear  in  the  x  or  y,- direction  and  t„  is 

X ,  y  d 

the  particle  arrival  time  in  hours.  The  first  term  on  the 
right  hand  side  of  Equation  III-38  accounts  for  toroidal 
growth  while  the  second  term  accounts  for  the  shear. 

The  toroidal  term  is  empirical  and  is  taken  from  WSEG- 
10  (4).  t*  is  equal  to  three,  hours  or  the  particle  arrival 
time,  whichever  is  less.  This  limit  on  t*  represents  the 
time  at  which  toroidal  motion  is  assumed  to  no  longer  have 
any  effect  on  the  cloud  spread  (12:7).  The  initial  cloud 
spread  oq  is  given  by  (25:51): 


a 

o 


1 .609EXI5 


.70  +  lnY/3 
m 


3.25 


4.0  +  ( InY  +  5.4)2 
m 


(III-39) 


where  Y  ^  is  the  yield  in  megatons  and  oq  is  in  kilometers. 
The  quantity  T  is  also  empirical  and  is  given  by  (25:51): 

i 

Tc  =  (Hc/5.0)  -  2.5(!!c/60.0)2  (111-40) 


where  T  is  in  hours  and  il  is  the  particle  injection 

« 

height  given  by  Equation  1 1 1  —  5  and  converted  to  kilofeet. 


The  shear  term  is  derived  in  Reference  26.  a  is  the 


cloud  spread  in  the  vertical  direction  and  is  given  by 
(25:5J ) : 


*z  =  (1. 609/5. 28)(.18Hc)  (III-41) 

where  a  z is  in  kilometers.  Sx>y  has  units  of  inverse 
hours . 

To  have  a  closed  dose  rate  contour.,  the  user-specified 
contour  value  must  be  located  twice  on  the  hotline.  One 
point  will  be  near  ground  zero  and  the  other  point  will  be 
located  somewhere  further  down  the  hotline.  These  points 
are  found  by  calculating  the  dose  rates  at  the  particle 
arrival  points  and  then  linearly  interpolating  between  the 
points.  The  dose  rate  at  ground  zero  is  assumed  to  be  zero 
for  calculational  purposes,  even  though  it  is  not  physical¬ 
ly  realistic  for  a  surface  burst. 

Dose  Rate  Contour  Mapping 

The  hotline  and  contour  points  are  stored  in  the  out¬ 
put  file  DOSERTE.  DOSERTE  is  the  input  file  to  program 
HYDMAP.  II Y DMA P  converts  the  coordinates  in  DOSERTE  from 
kilometers  to  degrees  longitude  and  degrees  latitude  and 
then  overlays  these  points  on  maps.  HYDMAP  automatically 
selects  the  range  in  longitude  and,  latitude  for  the  plots. 
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IV.  Evaluation  and  Results 


Overview 

The  theory  and  equations  presented  in  the  previous  sec¬ 
tion  are  very  similar  and  in  many  cases  identical  to  those 
used  in  another  variable  wind  smear  code,  REDRAM.  REDRAM 
was  written  by  Arthur  Hopkins,  and  was  the  first  fallout 

v'  I., 

code  to  use  the  spectral  coefficient  method. 

REDRAM  was  validated  against  data  from  the  Mt.  St. 
Helens  volcanic  eruption  and  from  a  small  conventional 
explosives  test  called  DIRECT  COURSE.  Mt.  St.  Helens  was 
taken  to  be  representative  for  a  debris  cloud  from  a  nucle¬ 
ar  explosion  in  the  megaton  range,  since  the  ash  cloud  rose 
to  a  height  of  22  kilometers.  In  the  DIRECT  COURSE  test, 
the  Defense  Nuclear  Agency  detonated  600  tons  of  high  explo¬ 
sives.  Thus,  the  data  was  representative  of  a  nuclear  explo¬ 
sion  in  the  kiloton  range.  Agreement  in  both  cases  between 
the  observed  and  REDRAM-calculated  hotline  locations  and 
particle  trajectories  was  excellent.  Hopkins  also  deter¬ 
mined  that  the  contribution  to  the  total  dose  rate  error 
from  the  spectral  coefficients  was  approximately  four  per¬ 
cent.  (15:57-58,112) 

The  accuracy  of  HYDRA  was  evaluated  by  comparing  HYDRA 
output  with  REDRAM  output  for  nine  test  cases.  The  test 
cases  are  summarized  in  Table  IV-1. 
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Table  IV-1 

Summary  of  Test  Cases  Used  In  Evaluating  HYDRA 


Shot 

# 

Yield 

(Kt) 

Location  of 
Deg.  Long. 

Ground  Zero 
Deg.  Lat. 

1 

500 

45.000 

-30.000 

2 

100 

354.000 

46.000 

3 

1000 

58.000 

75.000 

4 

100 

170.000 

-40.000 

5 

100 

253.525 

33.624 

6 

100 

165.230  ...  t.. 

11.350 

7 

1000 

37.600 

58.000 

8 

500 

90.000 

15.000 

9 

_ 

500 

300.000 

-60.000 

All  test  cases  in  Table  IV-1  used  the  following  data: 

1.  The  spectral  coefficients  were  derived  from 
"typical"  windg  for  January  14,  1978. 

2.  Median  particle  radius  =  0.204  microns. 
(DELFIC  default). 

3.  Logarithmic  slope  =  ln(4.0)  (DELFIC  default). 

4.  Fraction  of  the  particle  activity  located  in 
the  particle  volume  =  0.68  (from  DELFIC). 

5.  Particle  density  =  2600.0  kg/m^. 

6.  Fission  fraction  =  0.5. 

7.  Source  normalization  constant  =  6084.0 
roentgen-km Vhr-kt  (from  DELFIC). 

Plots  comparing  the  contour  maps  produced  by  each  code 
are  given  in  Figures  IV-1  and  IV-2  for  test  cases  2  and  3. 
Comparison  plots  for  the  remaining  cases  are  given  in  Appen 
dix  A.  Overall,  the  agreement  between  HYDRA  and  REDRAM  is 
excellent,  with  identical  hotlines  being  produced  in  all 
runs.  The  differences  in  X  and  Y  coordinates,  shear  terms, 
and  arrival  times  were  less  than  0.5%. 

In  all  runs,  the  calculation  of  the  point  where  a  con- 
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IV-1.  Test  Case  2:  HYDRA  vs.  REDRAM 


HYDRA 

REDRAM 


tour  crossed  the  hotline  produced  the  largest  differences 
between  the  two  models.  These  differences  are  summarized 
in  Table  IV-2,  The  distance  from  ground  zero  was  deter¬ 
mined  using  HYDRA  data.  The  error  listed  in  Table  IV-2  was 
calculated  by  dividing  column  2  by  column  3. 


Table  IV-2 

Largest  Differences  in  Jgst  Cases 


1 - — 

Shot 

» 

Largest  Distance 
(Km) 

Distance  from 

GZ  (Km) 

%  Error 

1 

37.426 

1211.109 

3.09 

2 

115.171 

2498.814 

4.61 

3 

40.807 

726.159 

5.62 

4 

64.355 

1494.749 

4.31 

5 

51.675 

1200.397 

4.30 

6 

16.398 

502.951 

3.26 

7 

58.387 

1773.009 

3.29 

8 

25.657 

747.419 

3.43 

9 

71.332 

1530.939 

4 . 66 

+ - + 


Differences  between  the  contour  maps  produced  by  the 
codes  arise  from  two  main  sources: 

1.  In  calculating  dose  rates,  REDRAM  models  the 
cloud  as  a  single  flat  pancake  cloud,  while 
HYDRA  models  the  cloud  as  a  set  of  pancake 
clouds,  one  cloud  for  each  particle  size. 

2.  Different  methods  are  used  to  calculate  g(t). 
Each  of  these  differences  is  discussed  in  one  of  the  follow¬ 
ing  sections. 

Single  Versus  Multiple  Pancake  Clouds 

i  ■  i 

At  very  high  and  very  low  yields,  HYDRA  and  REDRAM  do 
not  agree  well.  The  contours  produced  by  HYDRA  are  shorter 
and  fatter  than  those  produced  by  REDRAM.  Other  than  the 
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g(t)  function,  the  major  difference  between  the  two  codes 
lies  in  the  calculation  of  the  vertical  cloud  distribution 
azand  of  the  quantity  Tc  .  REDRAM  uses  the  single  cloud 
height  of  Equation  II-3  in  these  calculations,  while  HYDRA 
uses  a  different  cloud  height  for  each  particle  size  as 
determined  by  Equation  III-5. 

To  determine  the  effects  of  using,Lia  single  pancake 
cloud  (SPC)  as  opposed  to  multiple  pancake  clouds  (MPC),  a 
second  version  of  HYDRA,  called  HYDRAS,  was  created  that 
used  Equation  11-3  to  calculate  a  and  T,  .  Test  cases  for 
various  yields  at  tiie  location  of  test  case  6  were  run  on 
the  two  codes.  Results  are  d.isplayed  in  Appendix  B,  and 
show  that  HYDRA  and  HYDRAS  produce  nearly  identical  con¬ 
tours  for  yields  in  the  "critical"  range  of  10  kilotons  to 
2.0  megatons.  Outside  this  range,  however,  the  contours 
for  the  multiple  cloud  model  start  to  spread,  with  the 
spread  increasing  as  the  yield  gets  larger  or  smaller.  Th 
explanation  for  this  behavior  must  lie  in  the  calculation 
of  (Eqn  111-38),  because  it  is  the  only  calculation 

dependent  on  the  values  of  and  Tfi . 

<?2  as  a  function  of  arrival  time  is  plotted  in  Fig¬ 
ures  I V  —  3  through  1 V— 5  for  both  cloud  models  for  yields  of 
1,  500,  and  15000  kilotons,  respectively.  In  all  three 
figures,  the  value  of  a for  the  MPC  model  rises  rapidly 

from  a  low  value  to  some  limiting  value.  This  limiting 

* 

value  is  dependent  on  the  weapon  yield. 

For  the  yield  inside  the  critical  range,  the  limiting 
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Fig.  1V-3.  Vertical  Cloud  Distribution  vs.  Arrival  Time 
Fbr  a  1  Kiloton  Burst 
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Ficj.  IV- 5.  Vertical  Cloud  Distribution  vs.  Arrival  Tine 
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value  of  o  for  the  MPC  model  is  the  value  of  a  given  by 
the  SPC  model.  For  an  arrival  time  of  five  hours,  the  MPC 
o  ^  value  is  already  at  95  percent  of  the  SPC  value.  There¬ 
fore,  because  the  MPC  value  rises  so  rapidly  to  the  SPC 
value,  little  difference  is  seen  in  the  contours  for  yields 
in  the  critical  range,  especially  for  late  arrival  times. 
For  yields  outside  the  critical  r.gnge,  the  difference 

in  o  between  the  two  models  becomes  very  large  for  arrival 
z 

times  greater  than  five  hours.  The  ratio 'of  the  MPC 
value  to  the  SPC  value  for  the  1  kiloton  burst  ranged  from 
1.35  to  1.41,  while  that  of  the  15  megaton  burst  ranged 
from  1.43  to  1.54.  The  larger  o  value  results  in  much 

•  Z 

larger  deviations  in  the  horizontal  cloud  distrioution  (See 
Eqn  III-38)  for  the  MPC  model  as  compared  to  the  SPC  model. 
Thus,  at  any  given  point  on  the  hotline,  the  dose  rate  in 
the  MPC  model  will  be  less  than  that  in  the  SPC  model  (See 
Eqn  III-36),  resulting  in  a  shorter  contour  pattern  in  the 
MPC  model.  However,  since  the  deviations  in  the  cloud  dis¬ 
tribution  are  so  much  larger,  the  contour  patterns  in  the 
MPC  model  will  also  be  much  wider  than  those  for  the  SPC 
model.  Figure  B-l  in  Appendix  B  is  the  most  dramatic  exam¬ 
ple  of  the  differences  in  the  contour  patterns  of  the  two 
models. 

The  relationship  between  T  and  arrival  time  for  the 

l*  t 

two  models  is  very  similar  to  the  one  observed  between  oz 

* 

and  arrival  time.  However,  the  effects  of  large  differ¬ 
ences  in  Tc  are  small  compared  to  those  achieved  by  the 


differences  in  0 


z .  The  value  of  T  affects  only  the 
toroidal  term  of  Equation  III-38,  which  contributes  very 
little  to  the  horizontal  spread  of  the  cloud  for  arrival 
times  greater  than  three  hours. 

Calculation  of  g(t) 

Both  models  use  Equation  III-ll  to  calculate  g(t),  but 
employ  different  methods  to  determine  the  quantity  dr/dt. 
REDRAM  uses  Colarco’s  coefficients  to  determine  dr/dt, 
while  HYDRA  uses  Equation  III-12.  Plots  of  g(t)  versus 
arrival  time  are  given  for  runs  2  and  3  in  Figures  IV-6  and 
IV-7.  From  these  plots  it  is  seen  that  using  the  method  of 
Equation  I I I — 1 2  produces  a  lfiwer  value  for  g(t)  for  most 
arrival  times  than  do  Colarco’s  coefficients.  This  results 
in  lower  dose  rates  along  the  hotline  and  causes  the  con¬ 
tour  lines  to  lie  closer  to  the  hotline.  An  inspection  of 
Figures  IV-1  and  IV-2  and  the  plots  in  Appendix  A  shows 
that  the  HYDRA  contour  lines  lie  inside  the  REDRAM  contour 
lines  in  every  test  case. 

In  test  cases  2  and  3,  the  differences  between  the  con¬ 
tour  maps  were  quite  large  for  locations  far  from  the  burst 
site.  To  determine  if  the  g(t)  calculations  were  responsi¬ 
ble  for  these  differences,  a  third  version  of  HYDRA,  named 
HYDRAC  (C  for  Colarco),  was  created  that  used  Colarco's 
coefficients  to  calculate  g(t).  Test  'c'ases  2  and  3  were 
run  through  HYDRAC,  and  the  results  were  compared  to  the 
REDRAM  output.  These  results,  displayed  in  Figures  IV-8 
and  IV-9,  show  that  the  differences  between  the  two  models 
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have  been  almost  completely  eliminated.  Thus,  the  g(t) 
calculations  are  the  main  source  of  the  differences  observ¬ 


ed  in  the  contour  patterns  for  yields  ranging  from  10  kilo- 
tons  to  2  megatons.  Outside  this  range,  the  different 
cloud  models  used  in  the  two  codes  may  cause  large  differ¬ 
ences  in  the  contour  patterns. 

The  slight  differences  observed  j,jg(  Figure  IV-9  can  be 
eliminated  by  replacing  the  multiple  pancake  cloud  model  in 
HYDRAC  with  the  single  pancake  cloud  model.  This  was  done 
for  test  case  3,  and  the  contour  patterns  produced  by  the 
two  codts  were  indistinguishable  when  plotted. 

Changing  the  size  of  thq  "bounding"  particle  radii,  R2 
and  Rl,  in  Equation  III-12,  from  1.01  X  (radius)  and  0.99  X 
(radius)  to  1.005  X  (radius)  and  0.995  X  (radius),  respec¬ 
tively,  had  no  effect  on  the  value  of  dr/dt.  The  values  of 
dr/dt  were  the  same  out  to  four  or  five  decimal  places. 

In  addition  to  determining  dr/dt,  Colarco's  coeffi¬ 
cients  also  calculate  the  radius  of  the  particle  hitting 
the  ground  for  a  given  arrival  time.  Thus,  a  good  way  to 
check  the  accuracy  of  Colarco's  coefficients  is  to  compare 
the  Colarco  radii  with  the  radii  used  to  determine  the  hot¬ 
line.  Ideally,  the  two  sets  of  radii  should  be  identical. 
The  ratio  of  the  actual  particle  radius  to  the  Colarco  radl 
us  is  plotted  as  a  function  of  arrival  time  for  three  dif- 

i'  i 

ferent  yields  in  Figure  IV-10.  This  plot  shows  that,  for 

* 

yields  ranging  from  100  to  1000  kilotons,  Colarco's  coeffi¬ 
cients  almost  always  underestimate  the  particle  size  for 
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Fig.  IV- 10.  Actual- to- Calculated  Particle  Size  Ratio 
vs.  Arrival  Tine  for  Various  Burst  Sizes 


arrival  times  out  to  24  hours.  Colarco's  coefficients  are 
most  accurate  for  arrival  times  ranging  from  2.0  to  12.0 
hours  (accurate  to  within  2%),  with  differences  of  up  to  6% 
for  arrival  times  outside  this  range.  Thus,  all  other  para¬ 
meters  being  equal,  any  code  using  Colarco’s  coefficients 
will  be  less  accurate  than  HYDRA,  since  HYDRA  uses  the  actu¬ 
al  particle  sizes  to  calculate  dr/dt  £gther  than  empirical¬ 
ly  fit  coefficients. 


Optimum  Number  of  Particles 

The  optimum  number  of  particles  is  the  one  that  will 
reproduce  the  hotline  and  contours  with  a  fair  degree  of 
accuracy  while  using  as  little  computing  time  as  possible. 
This  optimum  number  was  determined  by  comparing  the  contour 
patterns  produced  by  10-,  15-,  and  20-particle  runs  with  25- 
particle  runs.  Comparisons  were  done  for  test  cases  2,3,4, 
5,7,  and  9.  Little  difference  was  observed  in  the  contour 
patterns  produced  by  the  15-,  20-,  and  25-particle  runs. 
However,  the  10-particle  runs  were  generally  at  consider¬ 
able  variance  with  the  other  runs  with  test  case  4  produc¬ 
ing  the  worst  agreement.  Thus,  it  can  be  concluded  that  at 
least  15  particles  are  needed  to  accurately  determine  the 
hotline  and  contour  patterns.  Contour  comparisons  for  test 
case  4  are  presented  in  Figures  IV-11  through  IV-13. 


Capabilities  and  Limitations 

HYDRA  and  HYDMAP  together  are  capable  of  producing  con¬ 


tour  patterns  overlayed  on  maps  anywhere  in  the  world. 
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Test  cases  1  through  9  are  plotted  with  maps  in  Figures  C-l 
through  C-9  in  Appendix  C. 

The  contour  patterns  and  hotlines  produced  by  HYDRA 
are  most  reliable  in  the  northern  hemisphere  for  latitudes 
ranging  from  10  to  80  degrees.  The  data  points  needed  to 
create  spectral  coefficients  are  sparse  at  the  equator,  and 
are  not  available  at  all  in  the  southern  hemisphere.  The 
test  cases  in  the  southern  hemisphere  used  wind  patterns 
that  are  mirror  images  of  the  ones  in  the  northern  hemi¬ 
sphere  .  (16) 

For  best  results,  the  weapon  yield  should  be  limited 
to  10  megatons  or  less.  Otherwise,  the  particle  injection 
height  is  above  the  highest  spectral  level,  and  winds  can 


no  longer  be  interpolated  between  the  spectral  levels. 
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V.  Summary >  Conclusions,  and  Recommendations 


user-specified  number  of  particle  size  groups  that  are  grav¬ 
ity-sorted  within  the  cloud.  All  particles  are  assumed  to 
be  spherical  in  shape,  and  are  transported  to  the  ground 
using  McI)onald-Da vies  fall  mechanics  in  a  discretely  layer¬ 
ed  atmosphere.  Horizontal  transport  of  the  particles  is 
accomplished  using  real  winds.  These  winds  are  represented 
as  sums  of  orthogonal  functions  in  spherical  coordinates, 
using  a  method  first  developed  at  the  National  Meteorologi¬ 
cal  Center. 


The  connection  of  the  particle  arrival  points  for  each 
particle  size  group  forms  the  hotline.  Contour  patterns 
about  the  hotline  are  determined  by  assuming  that  the  activ¬ 
ity  distribution  about  the  hotline  has  a  Gaussian  form, 
witli  the  peak  of  the  Gaussian  on  the  hotline. 

HYDRA  is  designed  for  local  fallout  prediction  (24 
hours  or  less).  It  is  capable  of  predicting  fallout  pat¬ 
terns  anywhere  in  the  world  except  for  the  poles.  However, 
the  best  results  will  be  achieved  by  limiting  the  burst 
locations  to  the  northern  hemisphere  from  10  to  80  degrees 
north  latitude.  No  wind  data  is  currently  available  for 
the  southern  hemisphere,  and  the  wind  vectors  at  the  equa- 


tor  are  not  reliable  due  to  the  scarcity  of  data  there.  At 
least  15  particles  are  needed  to  adequately  define  the  hot¬ 
line  and  contour  patterns. 

HYDRA  was  validated  against  REDRAM,  another  fallout 
code  using  spectral  methods  developed  by  Arthur  Hopkins 
(14,15).  Overall,  the  contour  patterns  produced  by  the  two 
codes  matched  very  well.  Differences  between  the  contour 
patterns  were  due  to  different  methods  used  to  calculate 
g(t)  and  to  the  different  ways  in  which  the  initial  stabil¬ 
ized  cloud  was  modeled. 

HYDRA  uses  a  finite  difference  method  to  calculate 
g(t),  while  REDRAM  uses  a  set  of  coefficients  derived  from 
particle  fall  mechanics  data.  The  method  used  in  HYDRA  is 
more  accurate.  The  particle  sizes  calculated  using  the 
coefficients  were  smaller  than  the  actual  particle  sizes, 
with  differences  ranging  up  to  six  percent.  The  g(t)  calcu 
lations  were  thfc  major  source  of  the  differences  observed 
in  Lite  contour  maps  of  the  two  codes. 

Different  methods  are  also  used  to  model  the  initial 
stabilized  cloud.  HYDRA  models  the  cloud  as  a  set  of  pan¬ 
cake  clouds,  one  for  each  particle  size  group,  with  each 
cloud  at  a  different  altitude.  R E D R A M  uses  one  pancake 
cloud  for  all  particle  size  groups.  Both  models  produce 
nearly  identical  results  for  weapon  yields  between  10  kilo- 
tons  and  2  megatons.  However,  outside  this  range  the  con¬ 
tour  patterns  produced  by  the  HYDRA  model  become  progres¬ 
sively  wider  and  shorter  compared  to  those  produced  by  the 


REDRAM  model.  The  HYDRA  model  is  more  realistic  in  that  it 
logically  predicts  that  particles  of  different  sizes  will 
be  injected  to  different  altitudes. 


Recommenda  tions 

1.  Determine  if  data  is  available  for  the  southern 
hemisphere.  If  this  data  is  available,  procedures 
need  to  be  devised  for  retrieving  it. 

2.  Modify  HYDRA  so  that  it  can  determine  dose  con¬ 
tours  as  well  as  dose  rate  contours.  The  needed 
data  item  is  already  present  in  the  code  to  deter¬ 
mine  t  h e  dose. 


3.  Hopkins  (15)  has  devised  a  method  to  compute  the 
dose  rate  at  arbitrary  distances  from  the  burst 
site.  This  method  should  be  reviewed  for  possible 
inclusion  in  HYDRA. 

4.  Develop  a  set  of  spectral  coefficients  for  those 
six  nuclear  tests  used  to  validate  DELF1C  and  oth¬ 
er  fallout  codes:  Johnnie  Boy,  Jangle-S,  Bravo, 
Small  Boy,  Koon,  and  Zuni.  This  data  would  be 
used  to  further  evaluate  HYDRA.  (15:119) 


*  mN  » 
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Dose  Rate  Contour  Comparisons 


HYDRA  vs.  REDRAM 


This  appendix  contains  seven  dose  rate  contour  compari¬ 
sons  between  HYDRA  and  REDRAM  for  Test  Cases  1  and  4-9 
described  in  Section  IV  of  this  report.  HYDRA  is  the  fall¬ 
out  code  developed  for  this  project.  REDRAM,  written  by 
Arthur  Hopkins  (14,15),  is  the  fallout  code  that  was  used 


to  evaluate  HYDRA. 
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Fig.  A- 2 .  Test  Case  4:  HYDRA  vs.  REDRAM 
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Pig.  A-6.  Test  Case  8:  HYDRA  vs.  REDRAM 


This  appendix  contains  nine  dose  rate  contour  compari¬ 
sons  between  HYDRA  and  HYDRAS  for  various  yields  at  the 
location  of  Test  Case  6.  HYDRA  uses  the  multiple  pancake 
cloud  method  to  model  the  radioactive  dust  cloud,  while 
HYDRAS  uses  a  single  pancake  cloud. 
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Fiq.  11-2.  HYDRA  vs.  HYDRAS  for  a  10  Mt  Burst 


Fig.  B-6,  HYDRA  vs.  HYDRAS  for  a  500  Kt  Burst 


160  E  162  E  161  E  166  E  168 


Fig.  B-8.  HYDRA  vs.  HYDRAS  for  a  10  Kt  Durst 


Fiy.  B-9.  HYDRA  vs.  HYDRAS  for  a  1  Kt  Burst 


Appendix  C 


HYDRA  Output  With  Maps 

This  appendix  contains  nine  HYDRA  dose  rate  contours 
overlayed  on  maps  for  Test  Cases  1  through  9  described  in 
Section  IV  of  tills  report.  The  yield  in  Test  Case  6  has 
been  changed  from  100  kilotons  to  15  megatons.  The  maps 
are  drawn  by  the  program  using  DISSPLA  software. 
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HYDRA  Output  for  Test  Case  1 


DOSE  RATE  CONTOURS 

YIELD  (KT)  -  1000.000  FISSION  FRACTION  -  .500 

RUED  (MICRONS)  -  .201  BETA  -  1.000 

GROUND  ZERO:  LONG (DEG.)  -  58.000  LOT (DEG.)  -  75.000 


Fig.  C-3.  HYDRA  Outixit  Cor  'Rest  Case  3 


DOSE  RATE  CONTOURS 

-  100.000  FISSION  FRACTION 

MICRONS)  -  .204  BETA  -  4.C 

ONGIOEG. )  -  170.000  LATIDEG. ) 


dose:  rate:  contours 

YIELD  (KT)  -  100.000  FISSION  ERRCTION  ■ 

RMED  (MICRONS)  -  .204  BETR  -  4.C 

GROUND  ZERO:  LONG  (DEG. )  -  -106.475  LRKDEG.) 


DOSE  RRTE  CONTOURS 

7 1  ELD  (KT)  -  1000.000  FISSION  FRRCTION  - 

RMED  (MICRONS)  -  .201  BETR  -  1.000 

GROUND  ZERO:  LONG (DEG.)  -  37.600  LRT ( DEG . )  - 

■ 


.5i 


c 
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UOSE  RATE  CONTOURS 

YIELD  (KT)  -  500.000  FISSION  FRACTION 

RUED  (MICRONS)  -  .201  BETA  -  4.C 

GROUND  ZERO:  10NG(0EG.)  -  90.000  LAT (DEG . ) 


o 


LEGEND 
°  -  HOTLINE 

«  -  10.000  ROENTGEN/HR 

*  -  I .000  ROENtGEN/HK 

•-  .100  ROENTGEN/HR 

*  -  .100  RUENfGEN/HR 

*  -  .010  ROENTGEN/HR 

»  -  .010  ROENTCEN/HR 


HYDRA  CXitDUt  for  Test  Case  3 


Appendix  D 


User's  Guide  For  Operating  HYDRA 


Introduction 

HYDRA  is  written  in  Fortran  V  and  is  designed  to  run 
on  the  Control  Data  Corporation  (CDC)  CYBER  750  computer. 
The  CYBER  currently  runs  version  2.4.2  of  the  Network  Opera 
ting  System  (NOS).  The  HYDRA  codes  also  use  the  AFIT 
procedural  file.  To  get  a  copy  of  this  file  for  your  CYBER 
account,  enter  the  following  set  of  commands: 

GET , PR0CFIL/UN=T8401 15 
SAVE , PROCFIL 
UPROC, PROCFIL 
BEGIN, PV,  ,  XXX 

here  XXX  is  any  three  characters  you  want  to  appear  on  the 
cover  page  of  your  printouts.  The  remainder  of  this  user's 
guide  is  divided  into  the  following  sections: 

1  .  Com  pi  Ling  II Y  DRA  . 

2.  Input  file  descriptions. 

3.  Running  HYDRA. 

4.  Output  file  descriptions. 

Compiling  HYDRA 

To  compile  HYDRA,  the  HYDRA  source  file  must  exist  as 
an  indirect  access  file  in  your  account.  The  source  file 
is  simply  called  HYDRA.  To  compile  HYDRA,  enter  the  follow 
ing  set  of  commands: 
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GET, HYDRA 
IMSL 

FTN5  ,  I=IIYDRA  ,  L=0,  B=IIYDRAB 
REPLACE,  IIYDRAB 
RETURN ,  * 


These  commands  will  create  a  compiled  version  of  HYDRA 
named  IIYDRAB.  The  IMSL  command  retrieves  the  IMSL  librar¬ 
ies  for  use  by  the  compiler,  and  allows  the  compiler  to 
recognize  the  IMSL  subroutines  used  in  HYDRA.  The  RP  com¬ 
mand  stores  IIYDRAB  as  an  indirect  access  file  in  your 
account . 


Input  File  Descriptions 


To  run  HYDRA,  the  following  files  are  needed  as  input: 


File  INITIAL 

File  INITIAL  must  exist  as  an  indirect  access 
fiLe  in  your  account.  It  consists  of  a  single  line 
containing  the  following  data  (in  order): 


1.  NP  -  The  number  of  particles  used  to  define 
the  hotline.  The  format  is  13.  At  least  15 
particles  should  be  used.  Up  to  99  particles 
may  be  specified,  but  it  is  recommended  that 
no  more  than  25  particles  be  used  to  keep  CPU 
time  as  low  as  possible. 

The  format  for  each  of  the  remaining  data  items  is 
F10.3. 

2.  Y LI)  -  The  weapon  yield  in  kilotons.  HYDRA 
will  be  most  accurate  for  yields  below  10 
megatons . 

3.  FF  -  The  fission  fraction.  FF  normally  has  a 
a  value  of  0.5. 

4.  N0RC0N  -  The  source  normalization  constant  in 
roentgen-km^/hr-kiloton .  One  value  for  this 


constant,  based  on  DELFIC  data,  is  6084.0. 

5.  RMED  -  The  median  particle  size  in  microns  of 
the  number-size  distribution.  The  DELFIC 
default  value  is  0.204. 

6.  BETA1  -  The  natural  log  of  BETA1  is  the  log 
slope  of  the  number-size  distribution.  The 
DELFIC  default  value  is  4.0. 

7.  FV  -  The  fraction  of  the  total  activity  that 
is  distributed  in  the  particle  volumes.  A 
typical  value  for  FV  is  0.68,  which  is  based 
on  DELFIC  data. 

8.  RII0F  -  The  particle  density  in  kg/m^.  This 
is  normally  set  to  2600  kg/m^. 

9.  DLONO  -  The  longitude  of  ground  zero  (GZ)  in 
degrees.  DLONO  ranges  from  0  to  360  degrees. 

10.  DLAT0  -  The  latitude  of  GZ  in  degrees.  DLAT0 
ranges  from  -85.0  to  85.0  degrees.  Using  a 
value  of  90  degrees  will  cause  HYDRA  to  blow 

up . 

11.  TEXIT  -  The  exposure  time  in  hours.  This 
data  item  is  not  used  at  present. 


File  CONTOUR 

File  CONTOUR  must  exist  as  an  indirect  access 
file  in  your  account.  This  file  contains  the  user- 
specified  contour  dose  rates  in  units  of  roentgcns/hr . 
Any  number  of  dose  rates  may  be  specified.  Each  dose 
rate  value  occupies  one  line.  The  format  for  each 
line  is  E10.3.  The  dose  rates  should  be  specified  in 
order  of  decreasing  value. 


File  LEVEL 

File  LEVEL  contains  the  altitude  levels  in  kilo¬ 
meters  for  the  12  spectral  layers.  Although  the  0.0 
altitude  level  is  NOT  a  spectral  layer,  it  must  be 


present  for  HYDRA  to  run.  There  is  one  altitude  per 
line.  The  format  for  each  line  is  F10.3.  This  file 
should  NEVER  be  altered  in  any  way.  File  LEVEL  must 
exist  in  your  account  as  an  indirect  access  file. 

File  TAPE72 

TAPE72  is  the  output  file  from  the  program 
AFGL4,  and  must  exist  as  a  direct  access  file  in  your 
account.  TAPE72  contains  12  sets  of  complex  spectral 
coefficients,  one  set  f..r  each  of  the  12  spectral 
layers.  There  are  three  complex  pairs  of  coefficients 
per  line,  and  the  format  of  each  line  is  6E13.7.  This 
file  should  NEVER  be  altered  in  any  way.  The  genera¬ 
tion  of  TAPE72  is  described  in  Appendix  F. 

Running  HYDRA 

Once  the  input  files  have  been  set  up,  you  are  ready 
to  run  HYDRA.  'First,  you  must  modify  the  AFIT  procedural 
file.  This  file  exists  in  your  account  as  an  indirect 
access  fiLe  named  PROCFIL.  Add  the  following  code  to 
PR0CF1L: 

.  PROC ,  G0I1YDRA  . 

RETURN,*. 

GET, HYDRAB, INITIAL, CONTOUR, LEVEL. 

ATTACH , SPEC0EF-TAPE7 2 . 

IMSL . 

HYDRAB. 

.  F.OR 
$ 

This  code  or  "proc"  should  be  added  immediately  after  line 
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number  50.  After  saving  PROCFIL,  type: 


BEGIN, PV, , XXX 

where  XXX  is  again  your  cover  page  designator.  This  will 
set  up  PROCFIL  to  run  the  "proc"  GOUYDRA. 

Type  the  command  GOUYDRA  to  execute  the  HYDRA  program. 
After  a  few  seconds,  HYDRA  will  flash  the  following  message 
on  the  screen: 

IF  YOU  WISH  TO  USE  FORMATTED  INPUT,  ENTER  "0"  TO  GO  ON. 
IF  YOU  WISH  TO  ENTER  DATA  FROM  THE  SCREEN,  ENTER  "1". 

Entering  "0"  next  to  the  program  prompt  (?)  tells  HYDRA  to 
use  the  data  from  file  INITIAL  in  determining  the  fallout 
pattern.  Entering  a  "1"  will  bypass  the  INITIAL  file  and 
allow  you  to  enter  data  from  the  screen  instead.  Enter  one 
value  for  each  program  prompt,  and  be  sure  to  use  decimals. 

When  HYDRA  has  completed  its  calculations,  the  word 
STOP  will  appear  on  the  screen,  and  the  NOS  system  prompt 
will  appear  (/).  HYDRA  takes  approximately  32  CPU  seconds 
to  do  a  25  particle  analysis. 

Output  File  Descriptions 

The  following  files  are  output  from  HYDRA: 

File  DOSE RTF 

File  DOSERTE  contains  hotline  and  contour  dose 
rate  values  and  geographic  coordinates.  There  is  one 
set  of  (X,Y)  coordinates  per  line.  The  format  for 
each  line  is  5E12.5,  except  for  the  first  line,  which 
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is  a  duplication  of  the  data  in  file  INITIAL  with  a 
few  additional  data  items.  Each  line  contains  the 
following  data  (in  order): 


1.  The  X  coordinate  in  kilometers  for  a  hotline 
or  contour  point. 

2.  The  Y  coordinate  in  kilometers  for  a  hotline 
or  contour  point. 

3.  Either  the  user-specified  contour  dose  rate 
or  the  calculated  dose  rate  on  the  hotline  in 
roentgens/hr . 

4.  Either  the  calculated  contour  dose  rate  (to 
compare  with  the  user-specified  dose  rate)  or 
the  hotline  dose  rate  in  roentgens/hr. 

5.  The  particle  radius  in  microns. 


This  file  should  not  be  modified  in  any  way.  It  is 
the  fiLe  used  by  the  program  IIYDMAP  to  map  the  hotline 
and  contour  lines.  To  save  this  file  for  future  use 
by  IIYDMAP,  enter  the  command: 


REPLACE, DOSERTE 


after  running  HYDRA. 

File  PARTICL 

File  PARTICL  contains  a  list  of  the  particle 
sizes  calculated  by  HYDRA  in  determining  the  hotline. 
The  particle  sizes  go  in  order  from  smallest  to  larg¬ 
est.  There  is  one  particle  size  per  line.  The  line 
format  is  F10.3.  This  file  is  not  used  by  any  other 
program. 
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Pile  D1AGN0S 


File  DIAGNOS  contains  diagnostic  data  for  each 
particle  size.  Some  of  these  variables  are  g(t),  the 
injection  height  of  the  particle,  horizontal  and  verti¬ 
cal  cloud  distributions,  and  shear  values.  All  data 
items  are  labeled  in  DIAGNOS,  and  the  units  are  speci¬ 
fied  in  the  source  code. 


Appendix  E 


User's  Guide  for  Operating  HYDMAP 


Introduction 

Like  HYDRA,  HYDMAP  is  written  in  Fortran  V  ..and  is 
designed  for  use  on  the  CDC  CYBER  750  computer.  The  AFIT 
procedural  file  is  also  used  by  HYDMAP  codes.  HYDMAP  takes 
the  file  DOSERTE  from  HYDRA  and  plots  the  coordinates  on 
world  maps  using  DISSPLA  software.  DOSERTE  must  exist  as 
an  indirect  access  file  on  your  account  for  HYDMAP  to  run. 
This  can  be  done  by  typing  the  command: 


REPLACE, DOSERTE 


after  running  HYDRA.  A  fiLe  description  of  file  DOSERTE  is 
given  in  Appendix  D  of  this  report. 


Running  HYDMAP 

To  run  HYDMAP,  the  AFIT  procedural  file  (named 
PR0CF1L)  must  have  the  following  "proc"  added  to  it: 


. PROC , MAPMAKE . 

RETURN,*. 

GET, HYDMAP, DOSERTE. 

FTN5,  l=IIYDMAP,  L=0 . 

ATTACH ,DISSPLA/UN=APPLIB. 

ATTACH, MAPDTA/UN=APPLIB. 

L  DS  ET , LI B= DISSPLA, PRES  ET=ZE RO . 

I.GO. 

ATTACH, CCLIB38/UN=LIBRARY. 

GET , CALPOP/UN=APPLI B . 

LIBRARY, DISSPLA, CCLIB38. 

CALPOP , 1038 , AF , XXX , INPUT , OUTPUT , TAPE99 , META . 
ROUTE, TAPE99 , DC=PU , UN=AF , ST=CSA . 

.  EOR 
.  * 
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where  XXX  is  your  initials.  This  "proc"  should  be  inserted 
after  the  "proc"  GOHYDRA.  MAPMAKE  automatically  retrieves 
the  files  and  libraries  necessary  to  run  IiYDMAP,  and  com¬ 
piles  and  executes  IIYDMAP. 

Once  the  MAPMAKE  "proc"  has  been  added  to  the  procedur¬ 
al  file,  type : 


BEGIN, PV, ,ZZZ 

where  ZZZ  is  the  three  letter  code  used  to  identify  your 
printed  output.  You  are  now  ready  to  execute  IIYDMAP. 

To  execute  IIYDMAP,  type  the  single  command: 

MAPMAKE 

This  executes  the  "proc"  listed  above.  After  a  few  sec¬ 
onds,  IIYDMAP  will  display  the  maximum  and  minimum  X  and  Y 
coordinates  (in  degrees)  it  has  calculated  for  the  axes  of 
the  plot,  along 'with  the  grid  spacing.  IiYDMAP  then  gives 
you  the  option  of  using  the  values  it  has  calculated  or 
choosing  your  own  values.  Encer  a  "0"  after  the  program 
prompt  (?)  to  use  the  calculated  values,  or  a  "1"  to  select 
your  own  values. 

If  you  choose  to  select  your  own  values,  IIYDMAP  will 
prompt  you  to  set  up  the  axes  values  and  grid  spacings,  one 
value  at  a  time.  Enter  only  one  value  after  each  program 
prompt,  and  be  sure  to  use  decimals.  The  values  for  the  X 
coordinates  must  be  specified  as  follows: 

1.  Western  Hemisphere:  0  to  -180  degrees. 
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2.  Eastern  Hemisphere:  0  to  180  degrees. 


Values  for  the  Y  coordinates  are  specified  in  a  similar 
manner  : 

1.  Northern  Hemisphere:  0  to  85  degrees. 

2.  Southern  Hemisphere:  0  to  -85  degrees. 

Because  IIYDMAP  uses  a  Mercator  projection  to  create  the 
maps,  a  value  of  90  or  -90  degrees  for  the  Y  direction  will 
cause  the  DISSPLA  software  to  fail.  For  any  plot,  the 
range  for  the  X  axis  should  at  least  equal  the  range  for 
the  Y  axis. 

Once  the  axes  and  grid  spacing  values  are  entered, 
IIYDMAP  will  display  on  the  screen  a  set  of  four-letter 
codes,  each  code  representing  the  map  of  a  section  of  the 
globe.  Choose  a  map  by  typing  its  code  next  to  the  program 
prompt.  BE  SUR'E  TO  ENCLOSE  THE  CODE  IN  SINGLE  QUOTES  OR 
THE  PROGRAM  WILL  FAIL. 

After  a  map  has  been  chosen  and  a  few  seconds  have 
elapsed,  some  DISSPLA  diagnostics  will  be  displayed  on  the 
screen.  A  few  more  seconds  will  pass,  and  the  program 
display  the  message: 

ENTER  POST-PROCESSOR  DIRECTIVES 

along  with  the  program  prompt.  Simply  hit  the  return  key 
here.  The  plot  is  then  sent  to  the  HARRIS  plotter  to  be 
printed.  This  printing  process  can  take  15  to  20  minutes. 
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Appendix  F 


User's  Guide  for  Creating  TAPE72 


TAPE72  is  the  input  file  to  HYDRA  which  contains  the 
spectral  coefficients.  To  create  TAPE72,  the  following 
files  are  needed:  ALLWIND,  MKTP71,  and  AFGL4.  Each  file  is 
described  in  one  of  the  following  sections. 


File  ALLWIND 


File  ALLWIND  contains  gridded  wind  data  for  the  12 
spectral  layers  for  the  following  month/days: 


1 . 

14 

January  1978 

2. 

17 

February  1977 

3. 

18 

March  1981 

4. 

11 

April  1981 

5. 

29 

May  1981 

6 . 

22 

June  1980 

7.  16  July  1981 

8.  28  August  1981 

9.  16  September  1983 

10.  24  October  1983 

11.  6  November  1982 

12.  26  December  1980 


This  file  should  never  be  altered.  File  ALLWIND,  after  it 
lias  been  renamed  TAPE40,  is  the  input  file  to  the  For  tan 
program  MKTP71.  ALLWIND  is  stored  in  the  Central  File 
System  of  the  Network  Operating  System  (NOS).  To  retrieve 
it  as  a  local  file  named  TAPE40,  enter  the  command: 


SREAD.TAPE40, ALLWIND, AC1=E840547 


File  MKTP71 

File  MKTP71  is  a  Fortrun  program  that  takes  the  data 
from  ALLWIND  as  Input  and  outputs  the  data  file  TAPE71. 
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TAPE71  contains  gridded  wind  data  for  the  user-selected 
month/day  in  a  form  which  can  be  read  by  the  Fortran  pro¬ 
gram  AFGL4.  MKTP71  was  written  by  Lt.  Jeff  Brown  of  Aero¬ 
nautical  Systems  Division,  Wright-Patterson  AFB,  Ohio. 

Fiie  AFGL4 

Fiie  AFGL4  is  a  Fortran  program  written  in  Fortran  IV. 
This  code  uses  the  gridded  wind  data  in  TAPE71  to  create 
the  spectral  coefficients  stored  in  TAPE72. 

Creating  "Proc"  MAKE72 

To  create  TAPE72,  first  insert  the  following  "proc"  in 
the  AFtT  procedural  file,  PROCFIL: 


. PROC , MAKE72 . 
GET,MKTP71 , AFGL4 . 
FTN5 ,  I=MK'fP71 ,  L=0 . 
LGO . 

RETURN , LGO , MKTP71 . 
REWIND , * . 
FTN,I=AFGL4,L=0. 
LGO. 

PURGE, TAPE72. 
REWIND, TAPE72. 

COPY , TAPE72 ,TP72A . 
RETURN, TAPE72. 
DEFINE, TAPE72. 
REWIND, TP72A. 

COPY ,TP72A ,TAPE72 . 
RETURN,*. 

.  EOR 
.  * 


Add  these  lines  after  "proc"  MAPMAKE.  Once  this  code  has 
been  added,  type  REPLACE , PROCFIL  to  store  the  new  version 
as  an  indirect  access  fiie.  Then,  type: 


BEGIN, PV, ,XXX 


where  XXX  is  your  three-letter  output  designator.  This 
command  updates  the  procedural  file  to  allow  you  to  run  the 
new  "proc"  MAKE72.  "Proc"  MAKE72  compiles  and  runs  MKTP71 
and  AFGL4. 

Creating  TAPE72 

To  create  TAPE72,  ensure  that  MKTP71  and  AFGL4  are 
present  as  indirect  access  files  on  your  account,  and  that 
TAPE40  (ALLWIND)  is  present  as  a  local  file.  Then,  type 
the  command: 


MAKE72 

In  sequential  order,  the  program  will  display  the  month, 
day,  and  year  for  a  set  of  gridded  winds.  To  select  a  set, 
type  a  "1"  after  the  program  prompt  (?).  Only  one  date  may 
be  selected  for  each  run  of  MAKE72.  Type  a  "2"  for  those 
sets  you  do  not  wish  to  use.  Then,  once  MAKE72  has  stopped 
running,  the  file  TAPE72  will  exist  as  a  direct  access  file 
in  your  account.  You  are  now  ready  to  run  HYDRA. 

If  a  TAPE72  file  already  exists  in  your  account,  it 
must  be  renamed  if  you  wish  to  preserve  it.  To  rename 
TAPE72,  type  the  following  set  of  commands  before  running 
MAKE72 : 

ATTACH, TAPE72 
DEFINE, XXXXXXX 
COPY .TAPE72 , XXXXXXX 
RETURN , * 

where  XXXXXXX  is  the  file  you  wish  to  store  the  data  in. 
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Appendix  G 


HYDRA  Source  Code 


This  appendix  contains  the  HYDRA  source  code.  HYDRA 
is  written  in  Fortran  V  and  uses  two  IMSL  subroutines, 
MDNOR  and  MDNRIS.  All  physical  variables  and  their  units 
have  been  defined  in  the  code. 


rjooooooooooooooooooooonooonoooooooooooooooooooooooooo 


HYDRA  -  A  VARIABLE  WIND  SMEAR  CODE 


HYDRA  (HOTLINE  YIELD-DEPENDENT  RESIDUAL  ACTIVITY)  FOLLOWS  THE  PATH 
OF  A  FALLOUT  PARTICLE  FROM  THE  HEIGHT  TO  WHICH  IT  IS  INJECTED  BY  A 
NUCLEAR  BURST  TO  ITS  FINAL  RESTING  PLACE  ON  THE  GROUND.  THE  DISTANCE 
BETWEEN  THE  POINT  OF  BURST  AND  THE  PARTICLE'S  LANDING  POINT  IS  DEPENDENT 
ON  THE  WINDS  PRESENT  IN  THE  AREA  AND  ON  THE  PARTICLE  SIZE.  GIVEN  THE 
SAME  WINDS,  SMALLER  PARTICLES  WILL  LAND  FARTHER  FROM  THE  BURST  POINT 
THAN  WILL  LARGER  PARTICLES.  THE  EFFECTS  OF  BOTH  AND  PARTICLE  SIZE  ARE 
ACCOUNTED  FOR  IN  THIS  CODE. 

THE  RATE  AT  WHICH  THE  PARTICLE  FALLS  IS  DEPENDENT  ON  PARTICLE  SIZE, 

AND  IS  CALCULATED  USING  THE  DA V IES-MACDONALD  EQUATIONS  FOR  SPHERICAL 
PARTICLES.  THIS  CODE  OPERATES  ON  THE  ASSUMPTION  THAT  ALL  PARTICLES  ARE 
SPHERICAL. 

VARIABLE  WINDS  ARE  TAKEN  INTO  ACCOUNT  BY  USING  THE  SPECTRAL  COEFFICIENT 
METHOD  DEVELOPED  BY  HOPKINS.  THE  ATMOSPHERE  IS  DIVIDED  IF  0  12  LAYERS. 

THE  WIND  DATA  IN  THIS  CODE  IS  SPLIT  INTO  X  AND  Y  DIRECTION  COMPONENTS. 

IT  IS  ASSUMED  THAT  THERE  IS  NO  WIND  COMPONENT  IN  THE  Z,  OR  VERTICAL 
DIRECTION. 

WRITTEN  BY  CAPT.  TONY  STRINES,  JR.  DURING  THE  PERIOD  OCT-DEC  1986. 

VARIABLE  DEFINITIONS 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 

13. 

14. 


15. 

16. 

17. 

18. 


AIRRHO  - 
AIRVIS  - 
BGTAl  - 

CC 

CCO 

CNUM  - 
COLRAD  - 
CU 

CUMAR  - 
CV 


DELZ  - 
DLATO  - 

DLONO  - 

DRATE  - 


DRDT  - 

DXKM  - 

DYKM  - 
EPS 


AIR  DENSITY  IN  KG/M3 . 

DYNAMIC  VISCOSITY  OF  THE  AIR  IN  KG/M-SEC. 

LOG ( BETA 1)  =  LOGARITHMIC  SLOPE  OF  THE  NUMBER-ST  IE 
DISTRIBUTION.  BETA1  =  4.0  FOR  DELFIC. 

SAME  AS  COLRAD. 

COLATITUDE  OF  GROUND  ZERO  IN  RADIANS. 

NUMBER  OF  DRATE' S  SPECIFIED. 

CO-LATITUDE  OF  THE  PARTICLE  POSITION  IN  RADIANS. 

COMPLEX  COEFFICIENT  FOR  THE  X  OR  LONGITUDINAL  WIND 
COMPONENT  IN  THE  SPECTRAL  METHOD. 

THE  CUMULATIVE  ACTIVITY-SIZE  DISTRIBUTION.  THE  ACTIVITY- 
SIZE  DISTRIBUTION  USED  IS  FREILING'S  APPROXIMATION. 

COMPLEX  COEFFICIENT  FOR  THE  Y  OR  LATITUDINAL  WIND 
COMPONENT  IN  THE  SPECTRAL  METHOD. 

THICKNESS  OF  AN  ATMOSPHERIC  LAYER  IN  KILOMETERS . 

LATITUDE  OF  GROUND  ZERO  IN  DEGREES.  THE  EQUATOR 
CORRESPONDS  TO  0  DEGREES  LATITUDE. 

LONGITUDE  OF  GROUND  ZERO  IN  DEGREES.  GREENWICH,  ENGLAND 
CORRESPONDS  TO  0  DEGREES  LONGITUDE. 

DESIRED  DOSE  RATES  IN  ROENTGENS/HR.  SPECIFIED  BY  USER  FOR 
THE  PURPOSE  OF  DETERMINING  THE  DOSE  RATE  CONTOURS.  MUST 
BE  LISTED  IN  THE  CONTOUR. DATA  FILE  IN  ORDER  FROM  THE 
HIGHEST  DOSE  RATE  TO  THE  LOWEST. 

TIME  RATE  OF  CHANGE  IN  THE  SIZE  OF  THE  PARTICLES  HITTING 
THE  GROUND  IN  MICRONS/HR. 

THE  DISTANCE  IN  KILOMETERS  TRAVELLED  BY  A  PARTICLE  AS  IT 
FALLS  THROUGH  ONE  ATMOSPHERIC  LAYER. 

SAME  AS  DXKM,  BUT  IN  THE  Y  DIRECTION. 

EPSILON  VALUE  USED  TO  CALCULATE  THE  ASSOCIATED  LEGENDRE 
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c 

c 

c 

c 

C  19. 
C 

C  20. 
C  21. 
C  22. 
C 

C  23. 
C  24. 
C 

C  25. 
C 

C  26. 
C  27. 
C  28. 
C  29. 
C  30. 
C 

C  3'. 
C 

C  32. 
C 

C  33. 
C  34. 
C  35. 
C 
C 

C  36. 
C 

C  37. 
C 

C  38. 
C  39. 
C 

C  40. 
C  41. 
C  42. 
C  43. 
C  44. 
C  45. 
C  46. 
C  47. 
C 
C 

C  48. 
C  49. 
C  50. 
C  51. 
C  52. 
C 

C  53. 


POLYNOMIALS  IN  THE  SPECTRAL  METHOD. 

EPS  =  SQRT((N*N  -  L*L)/(4*N*N  -  1)), 

WHERE  N  IS  THE  LONGITUDINAL  INDEX,  AND 
L  IS  THE  LATITUDINAL  INDEX. 

FF  -  FRACTION  OF  THE  WEAPON  YIELD  ARISING  FROM  FISSION 
REACTIONS 

FV  -  FRACTION  OF  THE  ACTIVITY  CONTAINED  IN  THE  PARTICLE  VOLUME. 

G  -  GRAVITATIONAL  CONSTANT  IN  M/SEC2. 

HCKF  -  HEIGHT  IN  KILOFEET  OF  THE  FALLOUT  CLOUD  IN  THE  PANCAKE 
CLOUD  APPROXIMATION  FROM  WSEG-10. 

HCKM  -  SAME  AS  HCKF  EXCEPT  IN  KILOMETERS. 

HL  -  LOWER  BOUND  LN  KILOMETERS  OF  THE  SPECTRAL  LAYER  INTO 
WHICH  THE  PARTICLE  FALLS. 

HU  -  UPPER  BOUND  OF  THE  SPECTRAL  LAYER  INTO  WHICH  THE  PARTICLE 
FALLS.  HU  IS  IN  KILOMETERS. 

JNJHGT  -  HEIGHT  OF  INJECTION  IN  KILOMETERS  FOR  A  PARTICLE. 

JCAP  -  RilOMBOIDAL  TRUNCATION  LIMIT  FOR  THE  SPECTRAL  METHOD. 

LCII  -  MIDPOINT  OF  AN  ATMOSPHERIC  LAYER  IN  KILOMETERS. 

I.VI.  -  SPECTRAL  LAYER  NUMBER,  STARTING  WITH  1. 

NKWVX  -  X  COMPONENT  OF  THE  VELOCITY  AT  (NEWX.NEWY)  OF  THE  NEW 
ALTITUDE  IN  METERS/SECOND. 

NEWVY  -  Y  COMPONENT  OF  THE  VELOCITY  AT  (NEWX.NEWY)  OF  THE  NEW 
ALTITUDE  IN  METERS/SECOND. 

NEWX  -  X  COORDINATE  OF  THE  PARTICLE  AFTER  IT  HAS  FALLEN  TO  ITS 
NEW  ALTITUDE. 


NEWY 

NEWZ 

Nl. 


NORCON  - 
Nl’ 

PRAI)  - 
PROB  - 

R24 

RADIUS  - 
RAD 

RADI.ON  - 
RADMIC  - 
R EARTH  - 


Y  COORDINATE  OF  THE  PARTICLE  AT  THE  NEW  ALTITUDE. 

BOUNDARIES  IN  KILOMETERS  FOR  THE  ATMOSPHERIC  LAYERS. 

NUMBER  OF  ATMOSPHERIC  LAYERS  INTO  WHICH  THE  ATMOSPHERE  IS 
DIVIDED.  THE  CENTER  OF  THE  HIGHEST  LAYER  IS  THE  INJECTION 
HEIGHT  OF  THE  PARTICLE. 

THE  SOURCE  NORMALIZATION  CONSTANT  FROM  DELF1C.  UNITS  ARE 
R0ENTGEN-KM2/IIR-K I LOTON . 

NUMBER  OF  PARTICLES  CHOSEN  TO  DETERMINE  THE  HOTLINE. 

GROUND  IN  24  HOURS. 

PARTICLE  RADIUS  IN  MICROMETERS. 

THE  VALUE  OF  THE  CUMULATIVE  ACTIVITY-SIZE  DISTRIBUTION 
FOR  THE  PARTICLE  SIZE  THAT  LANDS  AT  24  HOURS  (APPROXIMATE). 
RADIUS  IN  MICROMETERS  OF  THE  PARTICLE  THAT  WILL  FALL  TO  THE 
SAME  AS  RADMIC. 

SAME  AS  RADMIC. 

LONGITUDE  OF  THE  PARTICLE  POSITION  IN  RADIANS. 

PARTICLE  RADIUS  IN  MICROMETERS. 

RADIUS  OF  THE  EARTH  IN  MISTERS. 


RIIOF  -  DENSITY  OF  THE  FALLOUT  PARTICLES  IN  KG/M3. 

RMED  -  PARTICLE  RADIUS  IN  MICROMETERS  FOR  WHICH  HALF  OF  THE 

TOTAL  ACTIVITY  IS  CONTAINED  IN  PARTICLES  HAVING  A  RADIUS 
LESS  THAN  OR  EQUAL  TO  RMED.  RMED  =  0.204  FOR  DELF1C. 

RR  -  SAME  AS  RADEON. 

RRO  -  LONGITUDE  OF  GROUND  ZERO  IN  RADIANS. 

SI1R.X  -  SHEAR  IN  THE  X  OR  LONGITUDINAL  DIRECTION  IN  INVERSE  HOURS. 

SIIRY  -  SHEAR  IN  THE  Y  OR  LATITUDINAL  DIRECTION  IN  INVERSE  HOURS. 

SHX  -  SHEAR  IN  THE  X  OR  LONGITUDINAL  DIRECTION  IN  INVERSE 

SECONDS  FOR  A  GIVEN  ATMOSPHERIC  LAYER. 

SHY  -  SHEAR  IN  THE  Y  OR  LATITUDINAL  DIRECTION  IN  INVERSE  SECONDS 
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FOR  A  GIVEN  ATMOSPHERIC  LAYER. 

SIZE  -  SIZE  OF  THE  INTERVALS  ON  THAT  PORTION  OF  THE  CUMULATIVE 
LOGNORMAL  ACTIVITY-SIZE  DISTRIBUTION  WHERE  THE  PARTICLES 
ALL  LAND  INSIDE  OF  24  HOURS. 

SPCLEV  -  ALTITUDE  OF  A  SPECTRAL  LEVEL  IN  KILOMETERS,  CORRESPONDING 
TO  A  GIVEN  SET  OF  SPECTRAL  COEFFICIENTS. 

TARRJV  -  TIME  OF  ARRIVAL  IN  HOURS  FOR  A  PARTICLE  ON  THE  GROUND. 

TEX IT  -  TIME  OF  EXIT  IN  HOURS  FROM  A  CONTAMINATED  AREA. 

TFALL  -  TIME  LN  SECONDS  FOR  A  PARTICLE  TO  FALL  THROUGH  ONE  LAYER. 

TLA PS  -  ELAPSED  TIME  IN  SECONDS  SINCE  THE  BURST. 

VZ  -  VERTICAL  VELOCITY  OF  THE  PARTICLE  IN  M/SEC  AT  A  LAYER 
BOUNDARY. 

XO  -  X  POSITION  OF  GROUND  ZERO  IN  KILOMETERS.  X  =  0  AT  GREEN¬ 
WICH,  ENGLAND  AND  INCREASES  AS  ONE  GOES  FROM  EAST  TO  WEST. 

YO  -  Y  POSITION  OF  GROUND  ZERO  IN  KILOMETERS.  Y  =  0  AT  THE 
EQUATOR  AND  INCREASES  AS  ONE  MOVES  NORTH.  Y  DECREASES 
AS  ONE  MOVES  SOUTH  OF  THE  EQUATOR. 

YLI)  -  WEAPON  YIELD  IN  KILOTONS. 

^  a.  .j,  .J.  . i.  j.  j.  Oe  J.  .J.  ^  £  £  .j.  a,  £  ;J.  £  ...  .£££.>.£  jj.  ;  j.  Jj.  £  £  J.  a.  £  >}.  ;j_  £  >Jc  >);  5^  £  ^ 

PROGRAM  HYDRA 


c 

54. 

SIZE 

c 

c 

c 

c 

55. 

SPCLE 

c 

56. 

TARR1 

c 

57. 

TEX  IT 

c 

58. 

TFALL 

c 

59. 

TLA  PS 

c 

60. 

VZ 

c 

c 

61. 

XO 

c 

c 

62. 

YO 

c 

c 

c 

63. 

YLI) 

c 

X  < 

REAL  RADMIC,  NEWX(0:30),  NEWY(0:30),  NEWVX(0:30),  NEWVY(0:30), 

I  IIU,  IIL,  AVEVEL,  TEX  IT,  PRA1),  NEWZ,  SLOPE,  INTCPT,  SPCLEV(0:30) , 

1  G,  YLI),  RIIOF,  NEWVX1 ,  NEWVY1,  NEWVX2,  NEWVY2,  XO,  YO,  T24, 

1  TARRIV(3) ,  TFALL(0:30) ,  INJIIGT(3),  LYLI),  AIRRHO,  AIRVIS,  FV, 

1  VZ,  VZF(0:30) ,  RAD(3) ,  PI,  REARTH,  PROB,  HCKF,  HCKM,  LNYMEG,  R24, 
I  ALPHO,  BEI'A ,  CUMAR,  SIZE,  ALPII25,  RAl)IUS(lOO),  X,  IT,  HCKM2, 

1  SIIRX,  SIIRY ,  DEI.Z,  I.CII,  SIIX,  SHY,  RMEI),  BETA  1 ,  FF,  NORCON,  DRDT 

INTEGER  LVL,  JCAP,  ML,  NP,  1ER 

PARAMETER  (JCAP=30,NL=10) 

COMPLEX  CU( 1 2.JCAP+ I , JCAP+1 ) ,  CV( 1 2 , JCAP+1 , JCAP+l ) 

DOUBLE  PRECISION  COLRAD,  RADLON,  EPS( JCAP+1 , JCAP+1 ) ,  CCO,  RRO 

OPEN  (6,ITLE='II0TLINE'  ,ERR=100) 

C  CALCULATE  THE  EPSILON  COEFFICIENTS. 

CALL  EPSI.ON( JCAP, EPS) 

PI  =  3.141593 

REARTH  =  635G76G.0 

G  =  9.81 


C  READ  IN  THE  SPECTRAL  COEFFICIENT  DATA. 

OPEN  ( 5 , FI I.E= ' SPECOEF '  ,ERR=100) 

REWIND  5 
IX)  5  I  =  1,12 

READ  (5, ' (GEI3. 7,2X) ' )  ((CU(  I  ,N,I.) ,N=1  .JCAP+1 ) ,U1 , JCAP+1) 
READ  <  V  ( GE1 3 . 7 , 2X ) ' )  ( (CV( I , N , L) , N«1 , JCAP+ 1  ),!.=  !, JCAP+ 1 ) 
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5 


CONTINUE 
CLOSE  (5) 


C  READ  JN  THE  SPECTRAL  LEVEL  DATA. 

OPEN  (5, FJ LE=' LEVEL* ,ERR=100) 

DO  30  L  =  0,40 

READ  (5, ' (F10.3) ' ,END=40)  SPCLEV(I) 

30  CONTINUE 
40  CONTINUE 
CLOSE  (5) 

C  THE  USER  MAY  SPECIFY  EITHER  FORMATTED  OR  INTERACTIVE  INPUT. 

PRINT  V1F  YOU  WISH  TO  USE  FORMATTED  INPUT,  ENTER  "0”  TO  GO  ON 
PRINT  VJP  YOU  WISH  TO  ENTER  DATA  FROM  THE  SCREEN,  ENTER  "1".' 
PRINT*,'  ' 

READ  *,  1DATA 

IF  ( IDATA.NE.O)  GOTO  6 

OPEN  (5, FILE-' INITIAL' ,ERR=100) 

REWIND  5 

READ  (5, ' ( 13, I0F10.3) ' , ERR=100)  NP,  YLD,  FF,  NORCON, 

1  RMEI),  BETA1,  I'V,  RHOF,  DLONO,  DLATO,  TEXIT 
CLOSE  (5) 

GOTO  7 


PRINT  VENTER  THE  YIELD  IN  K1L0T0NS:' 

READ  *,  YLD 

PRINT  *,  'ENTER  THE  FISSION  FRACTION:' 

READ  *,  FF 

PRINT  *,  'ENTER  THE  SOURCE  NORMALIZATION  CONSTANT  (R-KM2/IIR-KT) 
READ  *,  NORCON 

PRINT  *,  'ENTER  THE  NUMBER  OF  PARTICLES:' 

READ  *,  NP 

PRINT  *,  'ENTER  THE  FALLOUT  PARTICLE  DENSITY  IN  KG/M3 : ' 

READ  *,  RHOF 

PRINT  *,  ' ENTER  THE  MEDIAN  PARTICLE  SIZE  IN  MICRONS:' 

READ  *.  RMED 

PRINT  *,  'ENTER  SLOPE  OF  THE  NUMBER-SIZE  DISTRIBUTION:' 

READ  *,  BETA  l 

PRINT  *,  'ENTER  THE  ACTIVITY- VOLUME  FRACTION:' 

READ  *,  FV 

PRINT  *,  ' ENTER  THE  LONGITUDE  IN  DEGREES  OF  GROUND  ZERO:' 


no 


HEAD  *,  DLONO 

PRINT  *,  'ENTER  THE  LATITUDE  IN  DEGREES  OF  GROUND  ZERO:' 
READ  *,  DLATO 

PRINT  *,  'ENTER  THE  EXIT  TIME  FROM  THE  AREA  IN  HOURS:' 
READ  *,  TEX  IT 

IF  (DLATO. GE. O.O)  THEN 

COLRAD  =  PI/2.0  -  DLAT0*PI/180.0 
ELSE 

COLRAD  =  -PI/2.0  -  DLAT0*P 1/180.0 
END  IF 

CCO  =  COLRAD 

RADLON  =  DLONO  *  PI  /  180.0 
RRO  =  RADLON 

YO  =  COLRAD  *  REARTII  /  1000.0 

XO  =  RADLON  *  REARTII  *  DSIN(COLRAD)  /  1000.0 

WRITE  (G,'(  14F10.3) ' )  YLD,  FV,  RMED,  BETA1 ,  RIIOF, 

1  DLONO,  DLATO,  RRO,  CCO,  XO,  YO,  TEX IT,  FF,  NORCON 


C  FIND  THE  DESIRED  NUMBER  OF  PARTICLES  THAT  WILL  ALL  STRIKE  TIIF,  GROUND 
C  INSIDE  OF  24  HOURS.  FIRST,  CALCULATE  THE  HEIGHT  OF  THE  PANCAKE  CLOUD 
C  USING  THE  WSEG-10  FORMULA. 

OPEN  (8 , F l  LE= ' PARTI  CL ' , ERR=100) 


LNYMEG  =  L0G( YLD/ 1000.0) 

IICKF  =  44.0  +  6. l*LNYMEG  -  .205*(LNYMEG+2.42)*ABS(LNYMEG+2.42) 
IICKM  =  ( 1 .009/5.28)  *  IICKF 
IICKM2  =  HCKM  /'2.0 


C 

C 

12 

15 


ESTIMATE  THE  S'ZE  OF  THE  PARTICLE  THAT  WILL  HIT  THE  GROUND  AT  EXACTLY 
24  HOURS  USING  MCDONALD-DAVIES  FALL  MECHANICS. 

CALL  STDATM( HCKM2 , A IRRIIO , A 1RVIS) 

R24  =5.0 


CALL  VPAI.UA IRRIIO, RIIOF, R24, A 1RV IS tVZ) 
T24  =  (HCKM* 1 000.0)  /  (VZ*3600.0) 

IF  (T24.LT. 24.0)  GOTO  15 
R24  =  R24  +  1.0 


GOTO  12 


R24  =  R24  -  1.0 


C  DEPENDENT  wlABLE  (.f  K.  THIS  CASE). 
A I. PI  10  =  I.OG(RMED) 

BETA  =  LOG(BETAl) 


C  THIS  IS  FREI LING'S  APPROXIMATION  OF  THE  FRACTIONATED  ACTIVITY-SIZE 
C  DISTRIBUTION. 

AI.PH25  =  A  I. PI  10  +  2.5  *  BETA— 2.0 
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TT  =  (l,OG(R24)  -  ALPII25)  /  BETA 
CALL  MDNOR  (‘FT,  PROB) 

SIZE  =  (1.0  -  PROB)  /  NP 

CUMAR  =  PROB  +  SLZE/2.0 

C  SUBROUTINE  MDNR1S  CALCULATES  THE  VALUE  OP  THE  DEPENEDENT  VARIABLE  (X 
C  IN  THIS  CASE)  OP  A  NORMAL  GAUSSIAN  DISTRIBUTION  GIVEN  THE  VALUE  OF 
C  THE  CUMULATIVE  PROBABILITY  OP  THE  GAUSSIAN  (CUMAR).  MDNRIS  IS  AN 
C  1MSL  SUBROUTINE. 

DO  10  I  =  1 ,NP 

CALI,  MDNRLS  ( CUMAR, X,1ER) 

RADIUS(l)  =  EX  PCX*  BETA  +  ALPII25) 

WRITE  (8, * (PJ0.3) ' )  RADIUS(I) 

CUMAR  =  CUMAR  +  SIZE 
10  CONTINUE 


CLOSE  (8) 


C  CALCULATE  THE  TIME  OP  ARRIVAL,  THE  FINAL  POSITION  ON  THE  GROUND,  AND 
C  THE  SHEAR  FORCES  FOR  EACH  PARTICLE  SIZE. 

OPEN  (5, PILE- 'PARTI CL' ,ERR=100) 

REWIND  5 

45  READ  ( 5 , ' ( PI 0 . 3) ' , END-90)  RADM1C 

RAD(l)  =  R A DM  1C  -  RADMIC*0.01 
RAI)(2)  =  RADMLC  +  RADM1C*0.01 
RAD(3)  =  RADMIC 

DO  75  K  -  1,3 

C  CALCULATE  THE  TIME  OP  ARRIVAL  FOR  EACH  PARTICLE. 

PRAD  =  RAI)(K) 

CALL  ARRTIM(  YLI),  PRAD.NL,  RIIOP,  DELZ,TA ,  PI  N  J I  IT ,  TFALL ) 

TARRIV(K)  -  TA 
INJIIGT(K)  =  PINJIIT 


C 


IP  (K.NE.3)  GOTO  75 


FIND  THE  NEW  X,  Y, 
NEWX(NUl)  =  0.0 
NEWY ( NL+ ! )  =  0.0 
SIIRX  =  0.0 
SHRY  =  0.0 
TUPS  =  0.0 
COLRAD  =  CCO 
RADLON  =  RRO 


AND  Z  COORDINATES  OF  THE  PARTICLE  AT  EACH  LEVEL. 


DO  70  I  =  NL,  ! .  -1 

‘FLAPS  =  ‘FLAPS  TFALL(l) 
NEW/.  =  I  *  DEL/. 

LCII  -  NEWZ  -  (0.5*DELZ) 
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c 

c 


81 


83 

C 


1 


1 


66 


CALCULATE  THE  AVERAGE  WIND  VELOCITY  ENCOUNTERED  BY  EACH 
PARTICLE  IN  EACH  ATMOSPHERIC  LAYER. 

CALL  LAYERS(  LCII ,  SPCLEV ,  LVL,  HU ,  IIL) 

IF  (LVL.EQ.l)  THEN 

CALL  UVCOMP( COLRAD , RADLON , LVL , EPS , CU , CV , NEWVX2, NEWVY2) 
NEWVX1  =0.0 
NEWVY1  =0.0 
GOTO  81 


END!  F 

IF  (LVL.EQ. 13)  THEN 
LVL  =  12 


CALL  UVCOMP(COLRAi) , RADLON , LVL , EPS , CU ,CV , NEWVX1 , NEWVY1 ) 
NEWVX(l)  =  NEWVX1 
NEWVY(l)  =  NEWVY1 
GOTO  83 
END  IF 


CALL  UVCOMP( COLRAD , RADLON , LVL, EPS , CU , CV , NEWVX2 , NEWVY2) 
LVL  =  LVL  -  1 

CALL  UVCOMP( COLRAD, RADLON , LVL, EPS , CU , CV , NEWVX1 , NEWV Y1 ) 


NEWVX(L)  =  (NEWVX2-NEWVX1)  *  ( LCII— IIL)  /  (IIU-IIL)  +  NEWVX1 
NEWVY(l)  =  (NEWVY2-NEWVY1)  *  (LCII-HL)  /  (HU-IIL)  +  NEWVY1 

IF  (I.EQ.NL)  GOTO  66 

CALCULATE  THE  ROOT-MEAN-SQUARE  SHEAR  IN  THE  X  AND  Y  DIRECTIONS. 
SNRX  =  SIIRX  +  ( ( NF.WVX ( 1+1 )  -  NEWVX(I))/(DELZ*1000.0))**2.0 

*  TFALL(H  l) 

SI  IX  =  SQRT(  SIIRX  /  TLA  PS) 

SIIRY  =  SIIRY  +  ((NEWVY(1+1)  -  NEWVY(i))/(DELZ*1000.0))**2.0 

*  TFAI.LU+1) 

SHY  =  SQRT(SIIRY  /  TLA  PS) 

CONTINUE 

DXKM  =  NEWVX(I)  *  TFALL(l)  /  1000.0 
NEWX(I)  =  NEWX( 1+1)  +  DXKM 

DYKM  =  NEWVY(l)  *  TFALL(l)  /  1000.0 
NEWY(I)  =  NEWY(Ul)  +  DYKM 


COLRAD  =  COLRAD  -  DYKM  *  1000.0  /  REARTII 
RADLON  =  RADLON  +  DXKM  *  1000.0 
1  /  (REARTII  *  1)81  N( COLRAD)) 

70  CONTINUE 


75  CONTINUE 


DR  DP  =  ADS(  (RAI)(I)  -  RAD(2))  /  (TARR1V(1)  -  TARR1V(2))  ) 
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El 

>  i 


SHRX  =  SQRT(  SIIRX*3600 . 0/TARRIV  ( 3) ) 

SI1RY  =  SQRT(SHRY*3600-0/TARRIV(3)) 

WRITE  (6,912)  RADMIC,  NEWX(l) ,  NliWY(l), 

1  COLRAD,  RADLON,  SIIRX,  S11RY,  DRDT,  INJIIGT(3),  TARRIV(3) 
912  FORMAT  (1  OK 12. 6) 

GOTO  45 

90  CLOSK  (5) 

Cl  JOSE  (6) 

C  CALCUI.ATK  THE  ACTIVITY  ALONG  THE  HOTLINE  AND  ALONG  THE  GAUSSIANS 
C  EXTENDING  FROM  THE  HOTLINE. 

CALL  ACTVTY 

100  STOP 
END 
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SUBROUTINE  ARRTIM 


c 

>Jc  $  $  $ : 

c 

c 

c 

1. 

ALT 

c 

c 

2. 

AVEVEI, 

c 

c 

c 

3. 

INTCPT 

c 

c 

4. 

PINJIIT 

c 

5. 

SLOPE 

c 

c 

c 

6. 

TA 

c 

7. 

TPALL 

c 

c 

8. 

VZF 

c 

V  X  A  A  X 

r  r  'r  r 

X  K  yK  A.  A  >r* 

v  V  'lv  X'  -»*  • 

&?* 


c 
c 

C  THIS  SUBROUTINE  CAI.CUl.ATES  THE  TIME  OF  ARRIVAL  IN  HOURS  FOR  A 
C  FALLOUT  PARTICLE  OF  ANY  SIZE. 


VARIABLE  DEFINITIONS 

-  THE  ALTITUDE  I.N  KILOMETERS  AT  AN  ATMOSPHERIC  LAYER 
BOUNDARY. 

-  AVERAGE  FALL  VELOCITY  IN  METERS/SEC  FOR  A  PARTICLE 
IN  A  GIVEN  ATMOSPHERIC  LAYER.  IT  IS  SIMPLY  THE 
AVERAGE  OF  THE  VELOCITIES  AT  THE  LAYER  BOUNDARIES. 

-  INTERCEPT  FOR  HOPKIN'S  EMPIRICAL  FITS.  UNITS  OF  MI 
(SEE  DTIC  //A  1 15— 514 ) . 


AS  A  FUNCTION  OF  PARTICLE  SIZE.  UNITS  OF 
METERS/MICROMETERS. 

-  TIME  OF  ARRIVAL  ON  THE  GROUND  IN  HOURS  FOR  A  PARTICLE. 

-  FALL  TIME  FOR  A  PARTICLE  IN  A  GIVEN  ATMOSPHERIC  LAYER 
IN  SECONDS. 

-  PARTICLE  FALL  VELOCITY  IN  METERS/SEC. 

s  ^ *  *  .*  *  -t-  *  *  *  *  -'n  --!<  *  *  *  *  -t-  *  *  *  *  *  *  *  *  *  #  *  *  *  *  *  *  *  *  *  *  *  $  *  #  *  *  *  ****#; 

SUBROUTINE  ARRTI M( YL1), I’RAD, NL, RIIOF, DELZ.TA , PINJIIT,  TFALL) 

REAL  I.YI.D,  INTCPT,  NEWZ,  VZF(0:30),  TFALL(0:30) 

INTEGER  NL 


C  CALCULATE  THE  INJECTION  HEIGHT  OF  THE  PARTICLE  USING  HOPKIN’S 
C  EMPIRICAL  FITS.  * 

I.YI.D  =  LOG(YLD) 

SLOPE  =  -KXP(  1.574  -  .OU97*LYLD  +  .03636*(LYLD**2.0) 

1  ~  .0041*(LYLD**3.0)  +  .0001965*(LYLD**4.0)) 

I NICPT  =  EX P( 7. 889  .34*I.YLD  +  .001226*(LYLD**2.0) 

1  -  . 005227* (LYl.D«'3.0)  +  .00041 7*(LYLI)**4.0)) 

PINJIIT  =  (SLOPE  *  2.0  *  PRAD  +  INTCPT)  /  1000.0 

C  HOPKIN’S  FITS  CAN  PRODUCE  NEGATIVE  INJECTION  HEIGHTS  FOR  LARGE 
C  PARTICLES  (>  300  MICRONS)  AND  LOW  YIELDS  (<  10  KT),  SO  THIS  SIMPLE 
C  FIX  WAS  ADDED  TO  INSURE  THAT  INJECTION  HEIGHTS  ARE  ALWAYS  POSITIVE. 

IF  (PINJHT.LT. 0.0)  PINJIIT  =  OLDINJ  *  OLDRAD  /  PRAD 
NEWZ  =  PINJIIT 

C  CALCULATE  THE  LAYER  THICKNESS.  THE  ATMOSPHERE  IS  DIVIDED  INTO  10 
C  LAYERS  WITH  THE  TOP  OF  THE  HIGHEST  LAYER  BEING  THE  PARTICLE  INJECTION 
C  HEIGHT. 

DEI.Z  =  NEWZ  /  (NL  -  0.5) 

NEWZ  =  NEWZ  +  (0.5  *  DELZ) 

C  CALCULATE  THE  AIR  DENSITY  AND  THE  DYNAMIC  VISCOSITY  AT  EACH  ALTITUDE 
C  LEVEL  (SUBROUTINE  STDATM) . 
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C  CALCULATE  THE  REYNOLDS  NUMBER  AND  THE  PARTICLE'S  FALL  VELOCITY  AT 
C  EACH  ALTITUDE  LEVEL  (SUBROUTINE  VFALL) . 

DO  50  r  =  NL,  0,  -1 
ALT  =  I  *  DELZ 

CALL  STDATM( ALT , AIRRHO , AIRVIS) 

CALL  VFALL(A1RRH0,RII0F,PRAD, AIRVIS, VZ) 

VZF(I)  =  VZ 
50  CONTINUE 


CALCULATE  THE  AMOUNT  OF  TIME  IT  TAKES  THE  PARTICLE  TO  FALL  THROUGH 
EACH  LEVEL. 


TA  =  0.0 

DO  60  1  =  NL,  1,  -1 

AVEVEL  =  (VZF(l)  +  VZF(I-l))  /  2.0 
TFALL(l)  =  1000.0  *  DELZ  /  AVEVEL 
TA  =  TA  +  TFALL(1 )/3600.0 
60  CONTINUE 


OLDINJ  =  PINJIIT 
OLDRAD  =  PRAD 
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C  SUBROUTINE  STDATM 

C 

C  TINS  SUBROUTINE  USES  THE  US  STANDARD  ATMOSPHERE  TO  CALCULATE  THE  AIR 
C  DENSITY  AND  DYNAMIC  VISCOSITY  AT  ALTITUDES  RANGING  FROM  0  TO  61,000 
C  METERS. 

0  $  Jic  'it  JSr'Jt-fcji:  j|- ft  ;|t  j|;  jfc  >};  ;J:  .J,t  £  £  ;|t  j^.  jJ:  jjc  .'it  jfc  >):  :i:  jit  :J:  :it  >it  Ji<  J|t  :|t  j[c  ;i:  &  ;it  J|t  J[t  Jit  J|t :}:  #  >\i  Jf:  >i;  if:  #  Jf;  >it  >}:  .'|t  J|t  Jit  Jf: 


c 

c 


VARIABLE  DEFINITIONS 


c 

1. 

AIRRHO  -  AIR  DENSITY 

IN  KG/M3  FOR  A  GIVEN  ALTITUDE. 

c 

2. 

AIRVIS  -  DYNAMIC  VISCOSITY  OF  AIR  IN  KG/M-S  FOR  A  GIVEN 

c 

3. 

ALT  -  ALTITUDE  IN 

METERS. 

c 

4. 

T  -  TEMPERATURE 

IN  DEG.  K  AT  ALT. 

c 

5. 

P  -  PRESSURE  IN 

PASCALS  AT  ALT. 

c 

6. 

PO  -  PRESSURE  IN 

PASCALS  AT  SEA  LEVEL. 
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SUBROUTINE  STDATM  (  NEWZ ,  A1RRII0 ,  AIR  VIS ) 

REAL  NEWZ,  AIRRHO,  AIRVIS,  P,  T,  PO,  ALT 

PO  =  1 .013E+05 

ALT  =  NEWZ  *  1000.0 

T  =  0.0 
P  =  0.0 

IF  (ALT.LT. 11000.0)  THEN 

T  =  288. IS  -  0.00654 5- ALT 
P  =  PO  *  (288.1 5/T )*#(-. 034 164/. 006545) 
GOTO  100 
END  IF 


IF  (ALT.GE. 11000.0  .AND.  ALT.LT. 20000.0)  THEN 
T  =  216.65 

P  =  22600.0  *  EXP( -.034164*( ALT-1 1000. 0)/216. 65) 
GOTO  100 
END  IF 

IF  (ALT.GE. 20000.0  .AND.  ALT.LT. 32000.0)  THEN 
T  =  216.65  +  .001  *  (ALT  -  20000.0) 

P  =  5528.0  *  (216. 65/T)**(. 034164/. 001) 

GOTO  100 
END  IF 

IF  (ALT.GE. 32000.0  .AND.  ALT.LT. 47000.0)  THEN 
T  =  228.65  +  .0028*( ALT-32000.0) 

P  =  888.8  *  ( 228. 65/T)**( . 034 164/. 0028) 

GOTO  100 
END  IF 

IF  (ALT.GE. 47000.0  .AND.  ALT.LT. 52000.0)  THEN 
T  =  270.65 


P  =  115.8  *  EXP(-. 034164  *  (ALT  -  47000. 0)/T) 
GOTO  100 
END  IF 

IF  ( ALT. GE. 52000.0  .AND.  ALT. LT. 61000.0)  THEN 
T  =  270.65  -  .002* (ALT  -  52000.0) 

P  =  62.21  *  (270. 65/T)**(-. 034164/. 002) 

END  1I? 

AiRRIlO  =  .003484  *  P  /  T 

A1RVIS  =  1.458E-06  *  (T**1.5)  /  (T  +  110.4) 
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SUBROUTINE  VEALL 


I’ll LS  SUBROUTINE  CM  Cl) BATES  THE  FALL  VELOCITY  OF  A  FALLOUT  PARTICLE 
AT  A  GIVEN  ALTITUDE  USING  DAVIES-MACDONALD  FALL  MECHANICS. 

£  #  £  # :‘f  sje  4; :};  $  £  ,-J;  *  ;1-.  jjs  ;j;  i{t  ;f;  jjt  £  j)t  s’;  jjc  $  $  #  $  ^  #  4: :):  ^  :{:  >[:  %  # }[:  :{:  :J:  ;J:  i}t  ^  i  #  >Y- 


VARIABLE  DEFINITIONS 


1.  RYNOLD  -  THE  REYNOLDS  NUMBER. 

2.  R2CD  -  THE  SQUARE  OF  THE  REYNOLDS  NUMBER  MULTIPLIED  BY  THE 

COEFFICIENT  OF  DRAG  OF  THE  PARTICLE. 


3.  RMET  -  RADIUS  OF  THE  PARTICLE  IN  METERS. 

4.  CF  -  EMPIRICAL  CORRECTION  FACTOR  FOR  THE  VELOCITY  OF  A 

SPHERICAL  PARTICLE  AT  11 J Gil  ALTITUDES.  TAKEN  FROM  DELFIC. 

5.  VZ  -  FALL  VELOCITY  OF  THE  PakTTCLE  IN  M/S. 


#  .f  .f.  $  -t.  ^  -h ^  ^  >n  $  £  %  $  * $  fc  $ & *  =!;  •■!:  $ # £  *  *  $  &  £  *  -'r  #  *  *  -t  *  -r  £ * #  ■'l''  *  *  *  *  *  *  *  *  ❖ 


SUBROUTINE  V  FALL  (A  J  RRIIO ,  RIIOF ,  RADIUS ,  A1RV1S ,  VZ) 

REAL  A  J  RRIIO,  RIIOF,  ALRVIS,  RYNOLD,  R2CD,  RMET,  G,  LIOREY,  RADIUS, 
J  LGR2CD 


G  =  9.81 

RMET  =  RADIUS  /  1.0E+06 


C 

C 


USE  MCDONALD-DAVJES  FALL  MECHANICS  TO  DETERMINE  THE  PARTICLE  FALL 
FALL  VELOCITY. 


R2CD  =  32.0  *  A 1  RRIIO  *  RIIOF  *  G  *  (RMET**3.0) 


/  (3.0*A1RVIS**2.0) 


IF  (R2CD.LE. 100.0)  THEN 

RYNOLD  =  R2CD/24.0  -  2.3363E-04*(lt2CD**2.0) 

1  +  2.0154E-06*(R2CD**3.0)  -  6.9105E-09*(K2CD**4.0) 

ELSE 

I.GR2CI)  -  L0G10(R2CD) 

LIOREY  =  -1.29536  +  .986*LGR2C!)  -  0.046677*(LGR2CD**2.0) 

1  +  .001 1235*(LGR2CD**3.0) 

RYNOLD  =  10.00**L10REY 
END  IF 

CF  =  1.0  +  I .  K)r)E-07/(RMET*AlRRlI0) 

VZ  =  RYNOLD  *  A1RV1S  *  CF  /  (2.0  *  AIRRIIO  *  RMET) 
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SUBROUTINE  EPSLON 


THIS  SUBROUTINE  CALCULATES  THE  EPSILON  COEFFICIENTS  FOR  COMPUTING 


THE  ASSOCIATED  LEGENDRE  POLYNOMIALS  WITH  THE  RECURSION  RELATION  IN 
NWS  30  P2A .  THE  CODE  IN  THIS  SUBROUTINE  IS  A  DUPLICATE  OF  THE  CODE 
FOUND  IN  IIOPKIN’S  MODEL. 


*  £  *  £  -t-  k  $  £  £  &  &  £  -t- 


X  X  A 
' i '  v  -i'  • 


A  A  .i.  .t  .1,  t,  A  A  .L  ,L  ,L  X  X  x  X  X  X  X  X  X 

^p»  rp  'I.  rf%  rp  .|\  rp  rp  .p  f|.  /p  .p  ^|S  rp  ^p  ^p  . 


VARIABLE  DEFINITIONS 


1. 

JCAP 

-  TRUNCATION  LIMIT. 

2. 

EPS 

-  EPSILON  COEFFICIENTS 

3. 

L 

-  LATITUDINAL  INDEX. 

4. 

N 

-  LONGITUDINAL  INDEX. 

**>! 


1  t  X  I.  A  A  A 


AtC  A  A  t 
T'  'P  -|'  V  'I'  •»*  ■ 


SUBROUTJ NE  EPSLON(JCAP.EPS) 

DOUBLE  PRECLS10N  EPS(l),  A 

JCAP1  =  JCAP  +  1 
JCAP2  =  JCAP  +  1 

DO  100  LI,  =  1.JCAP1 
L  =  LL  -  1 
JLE  =  L  *  JCAP2 

DO  110  1NDE  =  2.JCAP2 
N  =  L  +  INDE  -1 

A  =  (N-N  -  L*L)  /  (4.D+0  *  N*N  -  1.D+0) 
EPS (JLE+ INDE)  =  DSQRT(A) 

110  CONTINUE 

100  CONTINUE 

IK)  200  LL  =  1.JCAP1 

JLE  =  ( LL— 1 )  *  JCAP2 
EPS( JLE+1 )  =  0 . D+0 
200  CONTINUE 
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SUBROUTINE  UVCOMP 

THIS  SUBROUTINE  CALCULATES  THE  WIND  VECTOR  COMPONENTS  FROM  NWS 
SPECTRAL  COEFFICIENTS.  IT  DUPLICATES  A  SUBROUTINE  IN  HOPKIN’S 
CODE. 

THERE  ARE  JCAP+1  LATITUDINAL  INDICES  AND  JCAP+1  LONGITUDINAL  INDICES 
IN  THE  SPHERICAL  HARMONICS  SUMMATION. 

^  'p  'Jj  'c  V  »J*  »}»  V  »}»  »I«  V  ^  ^  ^  ^  *}*  *N  $  V  'I'  'p  '}*  ^  ^  »}»  ^  'p  'p  »}»  ^}»  ^  ^  ^  "Is  ^  ^ 


VARIABLE  DEFINITIONS 

-  COMPLEX  SPECTRAL  COEFFICIENTS  FOR  THE  X  (W-E)  WIND 
COMPONENTS. 

-  COMPLEX  SPECTRAL  COEFFICIENTS  FOR  THE  Y  (S-N)  WIND 
COMPONENTS. 

-  QUANTITY  EQUAL  TO  EXP(1L*LAMBDA) ,  WHERE  I  IS  THE 
J MAG J NARY  NUMBER,  L  IS  THE  LATITUDINAL  INDEX,  AND 
LAMBDA  IS  THE  LONGITUDE  IN  RADIANS. 

-  THE  X  OR  LONGITUDINAL  COMPONENT  OF  THE  WIND  VELOCITY  AT  A 
GIVEN  POSITION  IN  M/SEC. 

-  THE  Y  OR  LATITUDINAL  COMPONENT  OF  THE  WIND  VELOCITY  AT  A 
GIVEN  POSITION  IN  M/SEC. 

-  LATITUDINAL  INDEX. 

-  LONGITUDINAL  INDEX. 

8.  IJSUM  -  PSEUDO  WIND  COMPONENT  IN  THE  LONGITUDINAL  DIRECTION.  THE 
REAL  PART  OF  THE  SPECTRAL  EXPANSION. 

9.  VSUM  -  PSEUDO  WIND  COMPONENT  IN  THE  LATITUDINAL  DIRECTION.  THE 
REAL  PART  OF  THE  SPECTRAL  EXPANSION. 
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1. 

CU 

2. 

CV 

3. 

EIL 

4. 

US 

5. 

vs 

6. 

LL 

7. 

NN 

8. 

IJSUM 

9. 

VSUM 

:  5^  ❖  ❖  ■ 

f.  *.  -  x  v 

V  'i'  'i'  1'  'i”  ■ 

SUBROUT 1 NE  U VCOMP( CC , RR , LVL , EPS , CU , CV , US , VS) 
PARAMETER  (JCAP=30) 

COMPLEX  CU(I2, JCAP+1 , JCAP+1),  CV( 12, JCAP+1 , JCAP+1 ) , 
1  ML,  USUM,  VSUM 

DOUBLE  PRECISION  EPS(JCAP+1 , JCAP+1 ) ,  CC,  RR 

REAL  PLN(JCAP+1, JCAP+1) 

PI  =  3.141593 


JCAP1  =  JCAP  +  1 
JCAP2  =  JCAP  +  1 

CALL  PLN3(PLN,CC, JCAP, EPS) 

USUM  =0.0 
VSUM  =0.0 

DO  10  LL  =  l.JCAPl 
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L  =11,-1 

. . J5J  L  =..  CMl’l -X(  I.)CQS(UtKR)  .  DSIN.(L*RR.)  )..  -  . .  — 

DO  12  NN  =  1.JCAP2 
AA  =  2.0 

IF  (L.EQ.O)  AA  =  1.0 

tlSUM  =  USUM  +  AA  *  PLN(NN,LL)  *  CU(LVL,NN,LL)  *  EIL 
VSUM  =  VSUM  +  AA  *  PLN(NN,LL)  *  CV(LVL,NN,LL)  *  EIL 
12  CONTINUE 

10  CONTINUE 


IF  (AIJS(CC).LT.  l.D-10)  THEN 
PRINT*,'  ****  NORTH  POLE 
COTO  9999 


END  I  F 


C  CALCULATE  THE  WIND  VECTOR  COMPONENTS. 
US  =  REAL(USUM)  /  DSIN(CC) 

VS  =  REAL(VSUM)  /  DSIN(CC) 


9999  END 
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SUBROUTINE  PLN3 


THIS  SUBROUTINE  CALCULATES  THE  ASSOCIATED  LEGENDRE  POLYNOMIALS 
FOR  THE  WIND  COMPONENTS  USING  THE  RECURSION  RELATIONS  IN  BELOUSOV  AND 
NWS  30  P24 .  IT  IS  A  DUPLICATE  OF  A  SUBROUTINE  FOUND  IN  IIOPKIN'S  CODE. 


VARIABLE  DEFINITIONS 


1. 


2. 

3. 

Of  -V .*• 
'r'rv'i' 


PEN 


N 


-  ASSOCIATED  LEGENDRE  POLYNOMIAL.  THE  RECURSION  RELATION- 
USED  TO  CALCULATE  IT  IS: 

PLN(L,N+1)  =  (X*PLN(L,N)  -  EPS(L,N)*PLN(L,N-1)) 

/  EPS(L.N-H) 

WHERE  L  IS  THE  LATITUDINAL  INDEX, 

N  IS  THE  LONGITUDINAL  INDEX, 

PLN  IS  A  FUNCTION  OF  X, 

X  =  COS (LATITUDE  IN  RADIANS). 

-  LATITUDINAL  LNDEX. 

-  LONGITUDINAL  INDEX. 


>.  J.  ,V  .V  «■-  A  A  .V  A  A  A  A  A  .y  A  A  A  A  A  A  A  A 
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SUBROUTINE  PLN3(PLN,CC, JCAP.EPS) 

DOUBLE  PRECISION  CC,  EPS(l),  SJNLAT,  C0S2,  PROD,  A,  B,  FL, 
1  PI,  P2,  P3,  UNFLOW 


REAL  PLN(l) 

DATA  UNFLOW/O. 75D-73/ 

SINLAT  =  DCOS(CC) 

COS 2  =  I.D+O  -  SINLAT*SJNLAT 

PROD  =  l.D+O 

A  =  l.D+O 

B  =  O.DlO 

JCAP1  =  JCAP  +  I 

JCAP2  JCAP  +  1 

DO  300  LI.  =  I, JCAP  1 


FL  =  L 

JLK  .  I.  *  JCAP2 
IF  (L.EQ.O)  GOTO  400 
A  ■  A  +  2.D+0 
B  =  B  +  2.D+0 

C  FIX  ANY  UNDERFLOW  VALUES. 

IF  ( PROD. I.E. ONFLOW)  PROD  =  O.D+O 
PROD  =  PROD  *  COS 2  *  A  /  B 
400  CONTINUE 

PI  =  l)S()RT(0.riI)+0  *  PROD) 
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PI.N(  JI.R+I )  =  SNGL(PI) 

i>2 '=  usqrk i*(2: D+c/  •+  3:D+or-  s'Inlat  *  pi  *"  ' 
PLN(JU'+2)  =  SNGL(P2) 

DO  500  N  =  3, JCAP2 
I  .INDEX  =  JLE  +  N 

P3  =  (SINLAT*P2  -  EPS(LINDEX-1)*P1)  /  EPS(LINDEX) 
PI.N(L  INDEX)  =  SNGL(P3) 

PI  =  P2 
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SUBROUTINE  L, AYERS 


THIS  SUBROUTINE  LOCATES  THOSE  ATMOSPHERIC  LAYER  BOUNDARIES 
WHICH  BRACKET  THE  INJECTION  HEIGHT  OF  THE  PARTICLE. 

^  *  ;J;  -t-  *  *  >Jc  >J: ***:!:  :Jt  4:  sf:  *  *  *  *  #  $  >Je  ^  $  :Jc  >):  #  #  :{:  $  £  $  sje  #  # #  #  >!:  t-  $  sjt 


SUBROUTINE  LAYERS( LCH , SPCLEV , LVL, HU ,  1IL) 

REAL  LCH,  SPCLEV(0:30) ,  HU,  1IL 
INTEGER  LVL 

IF  (LCII.GT.SPCLEV(O)  .AND.  LCii.LE.SPCLEV(l))  THEN 
LVL  =  1 

HU  =  SPCLEV ( 1 ) 

111,  =  SPCLEV(O) 

GOTO  200 
END  IF 


LF  ( LCH. GT. SPCLEV ( 1)  .AND.  LCH.LE.SPCLEV(2))  THEN 
LVL  =  2 

IIU  =  SPCLEV(2) 

HL  =  SPCLEV(l) 

GOTO  200 
END  IF 


IF  (LCII.GT.SPCLEV(2)  .AND.  LCH.LE.SPCLEV(3) )  THEN 
LVL  =  3 

IIU  =  SPCLEV(3) 

III.  =  SPCLEV(2) 

GOTO  200 
END  IF 


IF  ( LCH .GT. SPCLKV( 3)  .AND.  LCH.LE.SPCLEV(A))  THEN 
LVL  =  A 

IIU  =  SPCI.EV(A) 

III,  =  SPCLEV(3) 

GOTO  200 
END  IF 


IF  (LCH . GT. SPCLEV ( 4 )  .AND.  LCH.LE.SPCLEV(5))  THEN 
LVL  =  3 
HU  =  SPCLEV(r.) 

HI.  =  SPCLEV (4) 

GOTO  200 
END  IF 

IF  ( LCH. CT. SPCLEV (3)  .AND.  LCH.LE.SPCLEV(6))  THEN 
LVL  =  () 

1111  =  SPCLEV((>) 

III.  ^  SI’Cl.EV(3) 

GOTO  200 


JJLj(I^U?kSPfJJJyf6}_  .AND.  LCH . LK . SPCLEV ( 7 ) )  THEN 


LVL  =  7 

IIU  =  SPCLEV(7) 
HL  =  SPCLEV(6) 
GOTO  200 
ENDIF 


IF  ( LCH . GT . SPCI.EV ( 7 )  .AND.  LC11.LE.SPCLEV(8))  THEN 
LVL  =  8 

HU  =  SPCLEV(8) 

HL  -  SPCLEV(7) 

GOTO  200 
ENDIF 

IF  ( LCH . GT . SPCLEV ( 8 )  .AND.  LCH.LE.SPCLEV(9))  THEN 
LVL  =  9 
HU  =  SPCLEV(9) 

HL  =  SPCLEV(8) 

GOTO  200 
ENDIF 

IF  ( LC1 ! . GT . S PC.LE V ( 9 )  .AND.  LCH.LE.SPCLEV(IO))  TIlfiN 
LVL  =  10 
HU  =  SPCLEV ( 10) 

HL  =  SPCLEV(9) 

GOTO  200 
ENDIF 

IF  (LCH.GT.SPCLBV(IO)  .AND.  LCH. I.E. SPCLEV ( i ! ))  THEN 
LVL  =  1 1 
HU  =  SPCLEV ( 1 1 ) 

HL  =  SPCLEV(IO) 

GOTO  200 
ENDIF 

IF  ( LCH. GT. SPCLEV ( 1 1 )  .AND.  LC1I.LE.SPCLEV(  12))  THEN 
LVL  =  12 
HU  =  SPCLEV( 12) 

HL  =  SPCLEV ( 1 1 ) 

GOTO  200 
ENDTF 


IF  ( I .CH . GT . SPCLEV (12))  THEN 
LVL  =  13 
HU  =  SPCI.EV(  12) 

IIL  =  SPCLEV  ( 12) 

ENDIF 


200  END 
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SUBROUTINE  ACTVTY 


THIS  SUBROUTINE  CALCULATES  THE  ACTIVITY  ON  THE  GROUND  GIVEN  THE 
HOTLINE  DATA.  IT  USES  A  SMEAR  EQUATION  WITH  A  GENERALIZED  F(X,Y) 
DISTRIBUTION  TO  CALCULATE  THE  DOSE  RATE  AND  DOSE. 

4  .■); if  if  )Ji  if >f  if  -I;  if  .'I".  if  -t- ;}:  £  if  :f  >f  f £ &  if £  $  if  if if  if  sj:  #  £  if $  $  .'Jc  -f  sjc  %  & 

VARIABLE  DEFINITIONS 


c 

c 

1. 

NORCON 

c 

2. 

ALPIIO 

c 

3. 

ALPII2 

c 

4. 

ALPII3 

c 

5. 

AR 

c 

6. 

BETA 

c 

c 

7. 

CXIIOT 

c 

c 

8. 

CXIIOT  l 

c 

9. 

CYIIOT 

c 

10. 

CYII0T1 

c 

11. 

DOSE 

c 

12. 

DRATE 

c 

c 

13. 

DRATEO 

c 

c 

14. 

DRDT 

c 

c 

15. 

DRMAX 

c 

c 

16. 

DSRTL 

c 

c 

17. 

DSRTR 

c 

c 

18. 

GTA 

c 

c 

c 

19. 

GXL 

c 

c 

20. 

CXR 

c 

21. 

GYI. 

c 

22. 

GYR 

c 

23. 

R EARTH 

c 

c 

24. 

SIIEARX 

c 

c 

25. 

SHEARY 

c 

c 

26. 

SICX 

c 

27. 

SIGY 

c 

28. 

SIGZ 

SOURCE  NORMALIZATION  CONSTANT  FROM  DELFIC.  UNITS  OF 
ROENTGEN-KM2/IIR-K 1 1  ,OTON . 

LN  OF  THE  MEDIAN  PARTICLE  RADIUS. 

LN  OF  THE  MEDIAN  PARTICLE  SURFACE  AREA. 

LN  OF  THE  MEDIAN  PARTICLE  VOLUME. 

ACTIVLTY-SIZE  DISTRIBUTION  IN  INVERSE  MICRONS. 

LOGARITHM LC  SLOPE  OF  THE  NUMBER-SIZE  DISTRIBUTION. 

X  POSITION  IN  KILOMETERS  FOR  A  USER  SPECIFIED  DOSE 
RATE  CONTOUR  VALUE  ON  THE  HOTLINE. 

BEGINNING  AND  END  X  LOCATION  IN  KILOMETERS  FOR  A  GIVEN 
DOSE  RATE  CONTOUR. 

SAME  AS  CXIIOT  BUT  FOR  Y  POSITION. 

SAME  AS  CXHOT1  BUT  FOR  THE  Y  DIRECTION. 

DOSE  AT  GX,  GY  IN  ROENTGENS. 

DOSE  RATE  CONTOUR  VALUE  IN  R/IIR  SPECIFIED  BY  USER. 

DOSE  RATE  ON  THE  HOTLINE  IN  R/HR  AT  A  PARTICLE  IMPACT 
POINT. 

RATE  OF  CHANCE  IN  SIZE  OF  PARTICLES  HITTING  THE  GROUND 
IN  MICRONS/IIOUR. 

MAXIMUM  DOSE  RATE  VALUE  ON  THE  HOTLINE  FOR  ALL  PARTICLE 
SIZES. 

CALCULATED  DOSE  RATE  VALUE  IN  R/HR  AT  CXL.GYL.  SHOULD  BE 
CLOSE  TO  DESIRED  "DRATK"  VALUE. 

CALCULATED  DOSE  RATE  VALUE  IN  R/HR  AT  CXR.GYR.  SHOULD  BE 
CLOSE  TO  DESIRED  "DRATK"  VALUE. 

FRACTIONAL  RATE  OF  ACTIVITY  ARRIVING  ON  THE  GROUND  AT 
TIME  T  IN  INVERSE  HOURS. 

X  POSITION  OF  A  POINT  OFF  THE  HOTLINE  IN  KILOMETERS,  ON 
THE  LEFT  SIDE  OF  THE  HOTLINE  AS  YOU  WALK  ALONG  THE 
HOTLINE  FROM  GROUND  ZERO. 

SAME  AS  CXI.  EXCEPT  THE  POINT  IS  LOCATED  ON  THE  RIGHT 
SIDE  OF  THE  HOTLINE. 

SAME  AS  GXL  EXCEPT  FOR  THE  Y  COORDINATE. 

SAME  AS  GXR  EXCEPT  FOR  THE  Y  COORDINATE. 

RADIUS  OF  THE  EARTH  IN  METERS. 

COMPONENT  OF  SICX  DUE  TO  WIND  SHEAR.  UNITS  OF  SQUARE 
KILOMETERS. 

COMPONENT  OF  S1CY  DUE  TO  WIND  SHEAR.  UNITS  OF  SQUARE 
KILOMETERS. 

STANDARD  DEVIATION  IN  THE  X  DIRECTION  OF  THE  CLOUD 
ACTIVITY  DISTRIBUTION.  UNITS  OF  KILOMETERS. 

SAME  AS  SICX  BUT  IN  THE  Y  DIRECTION. 

STANDARD  DEVIATION  OF  THE  VERTICAL  DISTRIBUTION  IN  THE 
CLOUD  IN  KILOMETERS. 


oooooooooooo 


29. 

SI  GO  - 

30. 

SLP 

31. 

TASK  - 

32. 

TOROID  - 

33. 

XHOT  - 

34. 

YHOT  - 

STANDARD  DEVIATION  OF  THE  INITIAL  HORIZONTAL  CLOUD 
DISTRIBUTION  IN  KILOMETERS. 

SLOPE  OF  THE  STRAIGHT  LINE  EXTENDING  FROM  GROUND  ZERO 
TO  A  PARTICLE  IMPACT  POINT  ON  THE  HOTLINE. 

TIME  IN  HOURS  DURING  WHICH  TOROIDAL  MOTION  IS  IMPORTANT 
IN  SPREADING  THE  CLOUD. 

COMPONENT  OF  S1GX  AND  SIGY  DUE  TO  TOROIDAL  MOTION  OF 
CLOUD.  UNLTS  OF  SQUARE  KILOMETERS. 

X  POSIT  LON  IN  KILOMETERS  ON  THE  HOTLINE  FOR  A  GIVEN 
PARTICLE. 

SAME  AS  XHOT,  EXCEPT  IN  Y  DIRECTION. 


^  aa  a  a  a  a  a  a  A  a  a  a  ^  a  a  a  a  A  v  a  Ax  a  a  A  •**  A  A  a  a  A  A  A  A  A  a  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A 


SUBROUTINE  ACTVTY 


REAL  YLD,  FV,  A,  B,  C,  BETA  1 , 

I  NORCON,  SLGX,  SLGY,  PI,  LOGRAD,  MFACT,  ALPHO,  RAI)LON(IOO), 

1  BETA,  ALPII2,  ALPII3,  DRATE(IOO) ,  XO,  YO,  COLRAD(IOO), 
l  FF,  TEX  IT,  LNYM,  IIC,  TC,  S1GO,  SIGZ,  GTA ,  AR,  REARTII, 

1  DRATEO(IOO),  GXR(IOO),  GYR(IOO),  DLONO ,  DLATO,  QUAN1 , 

1  RAI)MIC(IOO),  XIIOT(IOO),  YHOT(iOO),  VXIIOT(IOO),  VYHOT(IOO), 

1  DSRTR(IOO),  DSRTL(LOO),  SHRX(IOO),  SliRY(lOO),  DRDT(IOO), 

1  INJIIGT(IOO),  TARRIV(IOO),  QI(IOO),  Q2( 100) , 

1  Q3(  100) ,  LOC,  CXIIOT,  CYI10T,  GXL(IOO),  GYL(IOO) 

INTEGER  CNIIM,  NP,  M 

C  FI!.!*:  DIAGNOS  CONTAINS  DIAGNOSTIC  DATA  FOR  EACH  PARTICLE  SIZE,  AND  IS 
C  NOT  USED  AS  INPUT  TO  ANY  OTHER  PROGRAM. 

OPEN  (<), F I LE=' DIAGNOS'  ,ERR=JOO) 

WRITE  (9,073)  YLD,  DLONO,  DLATO 
673  FORMAT  ('YIELD  (KT)  =  '  ,P10.3,5X,  'DLONO  =  \F8.3,5X, 

1  'DLATO  =  ' ,F8.3,/) 

OPEN  ( 0 , FI  I ,E= ' DOSERTE ' , ERR= 100) 

open  (3,FILE='II0TLINH'  ,ERR=100) 

REWIND  3 

READ  (5, ’ ( 14FI0.3) ' )  YLD,  FV,  RUED,  BETA  1 ,  RHOF, 

1  DLONO,  DLATO,  RRO,  CCO,  XO,  YO,  TEX  IT,  FF,  NORCON 
WRITE  (6, ' ( 14F10.3) ' )  YLD,  FV,  RMED ,  BETA  1 ,  RIIOF, 

1  DLONO,  DLATO,  RRO,  CCO,  XO,  YO,  TEXIT,  FF,  NORCON 
DO  10  1  =  1,100 

READ  ( 3, ' ( 10E12.6) ' , END=20)  RADMlC(l),  XHOT(I),  YIIOT(I), 

1  COLRAD(l),  RADLON(I),  SIIRX(I), 

1  SllRY(l),  DKDT(i),  lNJHGT(I),  TARR1V(1) 

10  CONTINUE 

20  NP  =1-1 

CLOSE  (3) 


C  READ  IN  THE  USER-SPECIFIED  CONTOUR  DOSE  RATES 
OPEN  ( 3, FI l.K=’ CONTOUR' ,ERR=IOO) 
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30 

40 


REWIND  5 

DO  30  1  =  1,100 

READ  (5, 1 (E10.3) ' ,END=40)  DRATE(I) 

CONTJ NUE 
CNUM  =1-1 

PI  =  3.141593 
R EARTH  =  6356766.0 
ALPIIO  =  LOG(RMED) 

BETA  =  LOG ( BETA 1) 

ALPII2  =  ALPIIO  +  2.0  *  BETA**2.0 
ALPII3  =  ALPIIO  +  3.0  *  BETA** 2.0 
1.NYM  =  L0G(  YLI)/ 1000.0) 

SJGO  =  1.609  *  EXP(0.7  +  LNYM/3.0 
1  -  3. 25/(4. 0  +  (LNYM  +  5.4)**2.0)) 

F I N  IT Y  =  -l.OE+25 
ZERO  =0.0 

WRITE  (6, '(5E12.5)')  ZERO,  ZERO,  ZERO,  ZERO,  ZERO 


C  CALCULATE  THE  DOSE  RATE  ON  THE 


HOTLINE  AT  THE  PARTICLE  ARRIVAL  POINTS. 


DRMAX  =0.0 
DO  50  I  =  NP, 1 ,-l 

MI-ACT  =  1.0/  (SQRT(2.0*P1)  *  BETA  *  RADM1C(1)) 
LOGRAD  =  L0G(RADM1C(1)) 


AR  =  FV  *  MI-ACT  *  EXP(-0.5  *  ( (L0GRAD-ALPH3)/BETA)**2.0)  + 

1  ( 1 — FV )  *  MEACT  *  EXP (-0.5  *  ( (L0GRAD-ALP1I2)/BETA)**2.0) 


OTA  =  AR  *  l)RDT(l) 

WRITE  (9, mG6)  RAI)MIC(I) ,  GTA ,  l)RDT(I) 

666  FORMAT  ('RADMIC  =  ' , F 10.3.5X, 'GTA  =  ',E12.5,5X, 

1  'DRDT  =  \B12.5) 


C  CALCULATE  S1GY  AND  S1GX  FOR  THE  WINDS. 

SIGZ  =  0.18  *  INJIIGT(l) 

IF  (T'ARRI V ( I  ).LT.3.0)  THEN 
TASK  =  TARRIV(I) 

ELSE 

TASK  =3.0 
END  11- 

IIC  =  (5.28/1.609)  *  INJHGT(l) 

TC  =  0.2  *  IIC  -  2.5  *  (HC/60.0)**2.0 
SI GY  =  SQRT(S 100**2.0  *  (I  +  8.0*TASK/TC) 

1  +  (SIOZ*SIIRY(  1  )*TARRI  V(  l)/2.0)**2.0) 

SIGX  =  SQRT(S 100**2.0  *  ( l  +  8.0*TASK/TC) 

1  +  (SlGZ*SHRX(i)*TARRlV(I)/2.0)**2.0) 


TOROID  =  S 100**2.0  *  (1.0  +  8.0*TASK/TC) 
SIIEARX  =  (S1GZ*SHRX( 1 )*TARR1V( I )/2.0)**2.0 
SHEARY  =  (SIGZ*SHRY(l)*TARRlV(l)/2.0)**2.0 
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o  o  o  r: 


QI(1)  =  NORCON  *  YU)  *  FF  *  GTA  /  SQRT(2.0*PI) 

Q2(l)  =  1.0  /  SQRT((SIGY  *  XIIOT(I)  /  TARRIV(I))**2.0 
1  +  (SIGX  *  YIIOT(I)  /  TARRIV(I))**2.0) 

Q3(I)  =  -0.5  /  ((S1GY*XH0T(1))**2.0  +  (S1GX*YH0T(1))**2.0) 


DRATEO(l)  =  Q1(J)  *  Q2(l) 

C  WRITE  OUT  MORE  DIAGNOSTIC  DATA. 

WRITE  (9,667)  SIGZ.SIGO.TC 

667  FORMAT  (5X,'SIGZ  =  * ,F10.3,5X, ’SIGO  =  ’,F10.3,5X, 

1  1 TC  =  1 ,F10.3) 

WRITE  (9,670)  SIIRX(l) ,SIIRY(I) ,TARRIV(I) 

670  FORMAT  (5X,’SilRX  =  '  ,E12.5,5X,  'SIIRY  =  ',E12.5,5X, 

1  'TARRIV  =  ' , EI2.5) 

WRITE  (9,669)  TOROID, SUEARX, SMEARY 

669  FORMAT  (5X, 'TOROID  =  ' ,E12.5,5X, 'SIIEARX  =  ',E12.5,5X, 

I  'SMEARY  =  ',E12. 5) 

WRITE  (9,672)  INJMGT(I) 

672  FORMAT  ( 5X ,' PARTICLE  INJECTION  HEIGHT  =  \F10.3) 

WRITE  (9,668)  Ql(i) ,Q2(I) ,Q3(1) 

668  FORMAT  (5X,'Q1  =  ' ,E12.5,5X, 'Q2  =  ',E12.5,5X, 

1  ' Q3  =  ' ,E12.5) 

WRITE  (9,671)  S I GX , S LGY , DRATEO( I ) 

671  FORMAT  (5X,'SIGX  =  ' ,E12.5,5X, 'SIGY  =  \E12.5,5X, 

1  'DRATE  =  ' ,E12.5,/) 

IF  (DRATEO(I).GT.DRMAX)  THEN 
DRMAX  =  DRATEO(I) 

ENDI F 

WRITE  (6, ' (5E12.5) 1 )  XIIOT(I),  YilOT(I),  DRATEO(l) ,  DRATEO(I) , 
I  RADM1C(1) 

50  CONTINUE 


C 

C 

C 


THIS  LINE  SERVES  AS  A  MARKER  TO  THE  PLOTTING  PROGRAM  TO  INDICATE 


BREAKS  BETWEEN  HOTLINE  POINTS  AND  CONTOUR  POIN'IS 
DIFFERENT  SETS  OF  CONTOUR  POINTS. 

WRITE  (6, 1 (5E12.5) ’ )  ITNITY,  FI  NIT Y,  FINLTY, 


,  AND  ALSO  BETWEEN 
FINITY,  ITNITY 


DO  60  I  =  1 ,CNUM 


IF  (DRATE(l).GE. DRMAX)  GOTO  60 


THE  POINT  CXH0T1 , CYH0T1  IS  THE  STARTING  POINT  FOR  A  GIVEN  DOSE 
RATE  CONTOUR.  THIS  POINT  IS  FOUND  BY  LINEARLY  INTERPOLATING 
BETWEEN  GROUND  ZERO  AND  THE  HOTLINE  VALUE  CLOSEST  TO  GROUND 
ZERO.  THE  DOSE  RATE  AT  GROUND  ZERO  IS  TAKEN  TO  BE  ZERO. 


CXIIOTI  =  ( DRATE(  I  )/l)RATEO(NP) )  *  XHOT(NP) 
CYH0T1  =  ( l)RATE(  1  )/DRATEO(  NP)  )  *  YHOT(NP) 
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WRITE  (6, '  (5E12.5) ' )  CXII0T1 ,  CYIIOT1,  DRATE(T), 
1  DRATE(l) ,  ZERO 


C 

C 


FIND  THE  DOSE  RATE  ALONG  THE  GAUSSiANS  EXTENDING  FROM  THE 
ARRIVAL  POINTS  ALONG  THE  HOTLINE. 

DO  70  K  =  NP.1,-1 


1 

1 

1 


IF  (DRATEO(K).LE.DRATE(I))  GOTO  70 

QUAN1  =  LOG ( DR ATE (J )  /  (Q1(K)*Q2(K)))  /  Q3(K) 

IF  (XIIOT(K)  .NE.O.O)  SLP  =  YIIOT(K)  /  XIIOT(K) 

IF  ( YHOT(K) .LT.0.0  .AND.  XHOT(K) .EQ.O.O)  SLP  =  -1.0E+30 

IF  (YIIOT(K)  .GT.0.0  .AND.  XIIOT(K). EQ.O.O)  SLP  =  1.0E+30 

IF  (YIIOT(K). NE.O.O)  THEN 

A  =  Y1I0T(K)*YI10T(K)  +  2.0*X1I0T(K)*XII0T(K) 

+  XIIOT(  K)**4 . 0/YII0T(  K)**2.0 
B  =  -(2.0*XHOT(K)*YII0T(K)*YII0T(K)  +  4.0*XIIOT(K)**3.0 
+  2.0*XHOT(K)**5.0/YI10T(K)**2.0) 

C  =  2.0*XII0T(K)**4.0  +  XII0T(K)**2.0  *  Y1I0T(K)**2.0 
+  X!!0T(K)**6.0  /  YIIOT(K)**2.0  -  QUAN1 

END  IF 


IF  (SLP. LT.0.0)  THEN 

IF  (YHOT(K). GT.0.0)  THEN 
GXR(K)  =  (-B  +  SQRT(B*B 
GXL(K)  =  (-B  -  SQRT(B*B 

ELSE 

GXR(K)  =  (-B  -  SQRT(B*B 
*CXL(K)  =  (-B  +  SQRT(B*B 
END  IF 
END  IF 


4.0*A*C))  /  (2.0* A) 
4.0*A*C))  /  ( 2 . 0*A ) 

4.0*A*C))  /  (2.0*A) 
4.0*A*C))  /  (2.0* A) 


IF  (SLP. GT.0.0)  THEN 

IF  (YHOT(K). GT.0.0)  THEN 
GXR(K)  =  (-B  +  SQRT(B*B 
GXL(K)  =  (-B  -  SQRT(B*B 

ELSE 

GXR(K)  =  (-B  -  SQRT(B*B 
GXL(K)  =  (-B  +  SQRT(B*B 
END  IF 
ENDIF 


4.0*A*C))  /  (2.0*A) 
4.0*A*C))  /  (2.0*A) 

4.0*A*C))  /  (2.0*A) 
4.0*A*C))  /  ( 2 . 0*A ) 


IF  (SLP. HQ. 0.0)  THEN 
GXR(K)  =  XIIOT(K) 
GXL(K)  =  XHOT(K) 
ENDIF 


IF  (SLP. NE.O.O)  THEN 

GYR(K)  =  ( -XIIOT( K )/ Y!IOT( K ) )  *  GXR(K) 

1  +  YHOT(K)  +  XIIOT(K)*XIIOT(K)/YHOT(K) 
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GYL(K)  =  ( -X i IOT ( K ) / Y HOT ( K ) )  *  GXL(K) 

1  +  YIIOT(K)  +  XIIOT(  K  )-XHOT(  K  )  /  YIIOT(  K ) 

ELSE 

GYR(K)  =  -SQR'l'(QUANl)  /  XIIOT(K) 

GYL(K)  =  SQRT(QUANl)  /  XHOT(K) 

END1F 


DSRTR(K)  =  QL(K)  *  Q2(K)  *  EXP(  Q3(K) 

1  *  (GXR(K)^YIIOT(K)  -  GYR(K)*XH0T(K))**2.0  ) 

DSRTL(K)  =  Q1(K)  *  Q2(K)  *  EXP(  Q3(K) 

1  *  (GXL(K)*YIIOT(K)  -  GYL(K)*XH0T(K))**2.0  ) 


C 

C 

C 


WRITE  OUT  THE  DOSE  RATES  THAT  ARE  ON  THE  RIGHT  SIDE  OF 


OE  THE  HOTLINE  AS  YOU  WALK  FROM  GROUND  ZERO  ALONG  THE 
HOTLINE. 

WRITE  (6, ' (5E12.5) ' )  GXR(K) ,  GYR(K) ,  DRATE(I) , 

DSRTR(K) ,  RADMIC(K) 


70  CONTINUE 


C 

C 

C 


85 


80 


LOCATE.  THE  POINT  ALONG  THE  HOTLINE  WHERE  THE  DOSE  RATE  IS  THE 


USER-SPECIFIED  CONTOUR  DOSE  RATE.  THIS  IS  SIMPLY  A  LINEAR 
INTERPOLATION  BETWEEN  TWO  HOTLINE  DOSE  RATE  VALUES. 

DO  80  M  =  2,NP 

IF  (DRATE( I ) .GE.DRATE0(M-1 ) 

.AND.  DRATE( I ) . LT. DRATEO(M) )  GOTO  85 

GOTO  80 

LOG  =  (DRATE(l)  -  l)RATEO(M-l ))  /  (DRATEO(M)  -  DRATEO(M-l)) 
CXIIOT  =  LOG  *  (XHOT(M)  -  XIIOT(M-l))  +  XIIOT(M-l) 

CYIIOT  =  LOC  *  (YIIO'i'(M)  -  YIIOT(M-l))  +  YHOT(M-l) 

WRITE  (0, ' (5E12.5) ' )  CXIIOT,  CYIIOT,  DRATE(l) , 

DRATE(I) ,  ZERO 

OOTO  1\ 


CONTINUE 


C 

C 


IF  THE  USER  SPECIFIED  DOSE  RATE  IS  TOO  SMALL  TO  BE  ON  THE  HOTLINE, 
PUT  IN  A  MARKER  TO  INDICATE  THE  BREAK  IN  THE  CONTOUR. 


WRITE  ((), '  (5E12.5) ' )  UNITY,  FIN  IT  Y,  FINITY, 


FINITY,  FINITY 


C 

C 

75 


WRITE  OUT  THE  DUSK  RATES  THAT  ARE  ON  THE  LEFT  SIDE  OF  THE 
HOTLINE  AS  YOU  WALK  FROM  GROUND  ZERO  ALONG  THE  HOTLINE. 

DO  77  K  =  1 , NP 

IF  (DRATEO(K).LE.DRATE(I))  GOTO  77 

WRITE  («), '  (5E12.5) ' )  GXL(K) ,  GYL(K) ,  DRATE(l), 

DSRTL(K) ,  RADMiC(K) 


77 


CONTINUE 
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c 


I2NU  THE  DOSE  RATE  CONTOUR  AT  THE  POINT  WHERE  IT  STARTED. 


WRITE  (6, 1  (5E12.5) ' )  CX1IOT1,  CYIIOT1,  DRATE(l) , 

DRATE(T) ,  ZERO 

WRITE  (6, ' (5E12.5) ' )  F1N1TY,  F1N1TY,  F1NITY,  FINITY,  FINITY 


60  CONTINUE 


100  CLOSE  (6) 
CLOSE  (9) 
END 
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Appendix  11 


11YDMAP  Source  Code 

This  appendix  contains  the  IiYDMAP  source  code.  H  Y  DM  A  P 
is  written  in  Fortran  V  and  uses  DISSPLA  subroutines  to 
draw  the  contours  and  maps.  All  physical  variables  and 
their  units  have  been  defined  in  the  code. 
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•tW-  *  *  K*  m*  ***  * 


K*%Y 


s*.  yy  1  y 


X  ^  ^  ^  t  A  t  X  X  l>  X  O  X  X  X  *#  X  t  XX  x  »*#x  xxxxxx^*^x  XXXXXX  X  X  X  XXXXX  XXXX  XXXXXXX  vl^XXXXXXX 

<-|N  -|*  |  -'j'-  *£•  j  I  j'  jh  'p  p  f|t  ('■p'p  'pV'p^p'p'i''p*p'p'|''p'PY'p'|'T'ls'l''P'll'f,'l'T,i',pV,P'P'l''PT'P,i',rT'PVT'l''l',r'flTT'l1'r/l'TTrl*'i'T'il 

PROGRAM  HYDMAP 


c 

AND 

PLOTS  Til 

c 

THE 

HOTLINE 

c 

DATA 

AND  CON 

c 

X  A  X  ;i  X  X  X  X 

c 

WRITTEN  B 

c 

X,  X.  A  JU  A 

V  v  'r  f  r 

X  X  i,  X  X  X  X  A 

c 

c 

c 

1. 

AREA 

c 

2. 

BETA1 

c 

3. 

CCO 

c 

4. 

COAST 

c 

5. 

CX 

c 

6. 

CY 

c 

7. 

DEG  KM 

c 

8. 

DLATO 

c 

9. 

DLONO 

c 

10. 

DRATE 

c 

11. 

DX 

c 

12. 

DXMAX 

c 

c 

13. 

DX0R1G 

c 

c 

14. 

DXSTP 

c 

15. 

DY 

c 

16. 

DYMAX 

c 

c 

17. 

DYORIG 

c 

c 

18. 

DYSTP 

c 

19. 

FF 

c 

20. 

FV 

c 

21. 

NORCON 

c 

22. 

R EARTH 

c 

23. 

RIIOF 

c 

24. 

RMEI) 

c 

c 

25. 

RRO 

c 

26. 

TEX  IT 

c 

27. 

XMAX 

c 

c 

28. 

X0R1G 

c 

c 

29. 

XO 

c 

30. 

YLD 

c 

31. 

YMAX 

c 

c 

32. 

YORIG 

c 

c 

33. 

YU 

PROGRAM  "HYDRA" 
"HYDRA"  GIVES 
CODE  TAKES  THE 


THIS  PROGRAM  TAKES  THE  "DOSERTE"  FILE  CREATED  BY 
F,  DATA  OVER  A  USER  SPECIFIED  WORLD  MAP. 

AND  CONTOUR  POINTS  IN  KILOMETERS.  THIS 
AND  CONVERTS  IT  TO  DEGREES. 

-J- if &  *  £  -J:  -'i;  £  if  s|: if  s|:  :J:  if  f  if  if  if  jJ:  £  if  ;Jc  s':  .-J;  if  ;Jc  ^5  if  >}c  f  >):  >;<  if  if  if  if  $$$$ 

TIEN  BY  CAPT.  TONY  STRINES,  JR.  DURING  THE  PERIOD  NOV-DEC  1986. 

VARIABLE  DEFINITIONS 

THAT  AREA  OF  THE  WORLD  YOU  WISH  TO  MAP. 

LOG  SLOPE  OF  THE  ACTIVITY-SIZE  DISTRIBUTION. 

CO-LAT LTUDE  OF  GROUND  ZERO  IN  RADIANS. 

SPECIFIES  WHAT  COASTLINES  DISSPLA  IS  TO  DRAW. 

X  COORDINATE  IN  KILOMETERS  FOR  A  CONTOUR  POINT. 

SAME  AS  CX,  EXCEPT  FOR  THE  Y  COORDINATE. 

CONVERSION  FACTOR  TO  GO  FROM  KILOMETERS  TO  DEGREES. 
LATITUDE  IN  DEGREES  OF  GROUND  ZERO. 

LONGITUDE  IN  DEGREES  OF  GROUND  ZERO. 

USER  SPECIFIED  DOSE  RATE  CONTOUR  IN  R/IIR. 

X  COORDINATE  LN  DEGREES  LONG.  FOR  A  CONTOUR  POINT. 
MAXIMUM  X  COORDINATE  IN  DEGREES  LONG.  FOR  ALL  CONTOUR 
AND  HOTLINE  POINTS. 

MINIMUM  X  COORDINATE  IN 
AND  HOTLINE  POINTS. 

PLACE  A  GRID  LINE  EVERY 

Y  COORDINATE  IN  DEGREES 
MAXIMUM  Y  COORDINATE  IN 
AND  HOTLINE  POINTS. 

MINIMUM  Y  COORDINATE  LN 
AND -HOTLINE  POINTS. 

PLACE  A  GRID  LINE  EVERY  "DYSTP"  DEGREES  LAT. 

FISSION  FR ACT  LON  OF  WEAPON,  USUALLY  .5. 

FRACTION  OF  ACTIVITY  DISTRIBUTED  JN  PARTICLE  VOLUME. 
SOURCE  NORMALIZATION  CONSTANT  LN  R-KM2/IIR-K1LOTON. 

RADIUS  OF  THE  EARTH  IN  KILOMETERS. 

PARTICLE  DENSITY  JN  KG/M3. 

MEDIAN  PARTICLE  RADIUS  OF  THE  ACTIVITY-SIZE  DISTRIBUTION 
IN  MICRONS. 

LONGITUDINAL  POSITION  OF  GROUND  ZERO  IN  RADIANS. 

AMOUNT  OF  TIME  IN  HOURS  SPENT  IN  CONTAMINATED  AREA. 
MAXIMUM  X  COORDINATE  IN  KILOMETERS  FOR  ALL  CONTOUR  AND 
HOTLINE  PO«NTS. 

MINIMUM  X  COORDINATE  IN  KILOMETERS  FOR  ALL  CONTOUR  AND 
HOTLINE  POINTS. 

X  COORDINATE  JN  KILOMETERS  OF  GROUND  ZERO. 

YIELD  OF  THE  WEAPON  IN  KILOTONS. 

MAXIMUM  Y  COORDINATE  IN  KILOMETERS  FOR  ALL  CONTOUR  AND 
HOTLINE  POINTS. 

MINIMUM  Y  COORDINATE  IN  KILOMETERS  FOR  ALL  CONTOUR  AND 
HOTLINE  POINTS. 

Y  COORDINATE  IN  KILOMETERS  OF  GROUND  ZERO. 


DEGREES  LONG.  FOR  ALL  CONTOUR 

"DXSTP"  DEGREES  LONG. 

LAT.  FOR  A  CONTOUR  POINT. 
DEGREES  LAT.  FOR  ALL  CONTOUR 

DEGREES  LAT.  FOR  ALL  CONTOUR 
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PROGRAM  IIOTMAP 

REAL  YLI),  EV,  RMED,  BETA1,  RIIOF,  DLONO ,  DLATO ,  RRO,  CCO,  XO,  YO, 
1  TEX IT,  FF,  NORCON,  XMAX,  YMAX,  DX(IOOO),  DY(IOOO),  XORIG, 

1  YOR1G,  CX(IOOO),  CY(IOOO),  DRATE(IOOO).  PI,  REARTH,  DEGKM 


CHARACTER  AREA-4 ,COAST*4 


DIMENSION  LABEL(5),  IPKRAY(500),  LEGLAB(IO) 

PI  =  3.141593 

REARTH  =  6356.766 

DEGKM  =  (.5  *  PI  *  REARTH)  /  90.0 

OPEN  ( 9 , F l LE= ' CON  ERR ' , ER R= 100) 

OPEN  ( 5 , FILE- ' DOSERTE ' , ERR-100) 

REWIND  5 

READ  (5, '  ( 14F10.3) 1 )  YLD,  FV,  RMED,  BETA  1 ,  RIIOF, 

1  DLONO,  DLATO,  RRO,  CCO,  XO,  YO,  TEXIT,  FF,  NORCON 

IF  (DLONO. GT. 180.)  THEN 
DLONO  =  DLONO  -  360.0 
END  IF 


XMAX  =0.0 
YMAX  =  0.0 
XORIG  =0.0 
Y0R1G  =  0.0 

C  LOCATE  THE  MAXIMUM  X  AND  Y  VALUES  AND  THE  MINIMUM  X  AND  Y  VALUES 
C  (XORIG,  YORIG) . 

DO  5  1  =  1,1000 

READ  (5, '(3E12.5) ' ,END=10)  CX(I),  CY(I),  DRATE(I) 

IF  (CX( I ) .EQ.-l .OL+25)  GOTO  5 

IF  (CX(I).GT.XMAX)  THEN 
XMAX  =  CX(l) 

XMAXY  =  CY(l) 

END  IF 


IF  (CX(I).I.T.. XORIG)  THEN 
XORIG  =  CX(I) 

X0R1GY  =  CY( I ) 

END1F 


IF  (CY(I).GT.YMAX)  THEN 
YMAX  =  CY(l) 

END  IF 


IF  (CY(l).LT. YORIG)  THEN 
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YORLG  =  CY(L) 
ENDIF 


5  CONTINUE 

10  REWIND  5 

DXMAXY  =  DLATO  +  XMAXY  /  DEGKM 

DXMAX  =  DLONO  +  XMAX  /  (1)EGKM*C0S(DXMAXY*PI/180. )) 

DXORY  =  DLATO  +  X0R1GY  /  DEGKM 

DX0R1G  =  DLONO  +  X0R1G  /  (DEGKM*C0S(DX0RY*PI/180. )) 

DYMAX  =  DLATO  +  YMAX  /  DEGKM 
DY0R1G  =  DLATO  +  Y0R1G  /  DEGKM 

IP  ( DXMAX. GE. 0.0)  THEN 

DXMAX  =  N1NT(A13S(DXMAX)/10.  +  .5)  *  10. 

IP  (DXMAX. GT. 180.)  THEN 
DXMAX  =  DXMAX  -  360. 

END  IF 

ELSE 

DXMAX  =  NINT( ADS (DXMAX)/ 10.  -  .5)  *  10. 

DXMAX  =r  -DXMAX 
END1P 

IP  (DYMAX. GE. 0.0)  THEN 

DYMAX  =  N1NT( ADS (DYMAX)/ 10.  +  .5)  *  10. 

IP  ( DYMAX. GT. 90.)  THEN 
DYMAX  =  180.  -  DYMAX 
ENDJLF 

ELSE 

DYMAX  =  N INT( ADS( DYMAX )/ 10.  -  .5)  *  10. 

DYMAX  =  -DYMAX 
END1P 


IP  (DX0R1G.GE.0.0)  THEN 

DXORJG  =  N 1 NT ( A DS ( DXOR 1G ) / 10 .  -  .5)  *  10. 

ELSE 

DXOR IG  =  N1NT(ADS(DX0RIG)/10.  +  .5)  *  10. 
DX0R1G  =  -DX0R1G 
IP  (DXORIG.LT. -180.0)  THEN 
DXOR IG  =  360.0  +  DX0R1G 
END  IF 
END1F 


IP  (DY0R1G.GE.0.0)  THEN 

DYORIG  =  NINT(ADS(DY0R1G)/10.  -  .5)  *  10. 

ELSE 

DYORIG  =  N INT(ADS(DY0R1G)/10.  +  .5)  *  10. 
DYORIG  =  -DYORIG 
IF  (DYORIG. LT. -90.)  THEN 
DYORIG  =  180.0  +  DYORIG 
ENDIF 
END1F 
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DXSTP  =  ABS( DXMAX  -  DXORIG)  /  5. 

DYSTP  =  ABS(DYMAX  -  DYORIG)  /  5. 

PRINT  *, ’MAX  X  (DKG.  LONGITUDE)  =  DXMAX 
PRINT  VMIN  X  (DEG.  LONGITUDE)  =  '.DXORIG 
PRINT  *, 'GRID  LINES  EVERY  '.DXSTP,'  DEG.  LONGITUDE' 
PRINT-,'  ' 

PRINT  *, 'MAX  Y  (DEG.  LATITUDE)  =  ' , DYMAX 

PRINT  *,'M1N  Y  (DEG.  LATITUDE)  =  '.DYORIG 

PRINT  *, 'GRID  LINES  EVERY  ' , DYSTP, '  DEG.  LATITUDE' 

PRINT*,'  '' 

PRINT  *, ' IF  THESE  VALUES  ARE  ACCEPTABLE,  ENTER  0  TO  GO  ON.: 
PRINT  VIE  YOU  WANT  TO  PICK  YOUR  OWN  VALUES,  ENTER  1.' 


READ  *,  NOYES 
IF  (NOYES. EQ.O)  THEN 
GOTO  22 
ELSE 

PRINT  VENTER  MAX  X  (DEG.  LONG’.):' 

READ  *,  DXMAX 

PRINT  *, 'ENTER  MIN  X  (DEG.  LONG.):' 

READ  *,  DXORIG 

PRINT  *, 'GRID  LINE  EVERY  (DEG.  LONG.):' 
READ  *,  DXSTP 
PRINT*,'  ' 

PRINT  *, 'ENTER  MAX  Y  (DEG.  LAT.):' 

READ  *,  DYMAX 

PRINT  *, 'ENTER  MIN  Y  (DEG.  LAT.):' 

READ  *,  DYORIG 

PRINT  * , '  GR 1 1)  LINE  EVERY  (DEG.  LAT.):' 
READ  *,  DYSTP 
BND1F 


CONTINUE 


PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 


❖ 

1 

»•» 


ENTER  THE  MAP  YOU  WISH  TO  USE: ' 

"PAFR"  FOR  AFRICA,' 

"PAS l"  FOR  ASIA,' 

"PAUS"  FOR  AUSTRALIA,' 

"PFUR"  FOR  EUROPE,' 

"PNOR"  FOR  NORTH  AMERICA,' 

"POL I"  FOR  THE  ENTIRE  WORLD,' 

"PSOU"  FOR  SOUTH  AMERICA,' 

OR  "USAM"  FOR  UNITED  STATES  BOUNDARIES.' 

I 

BE  SURE  TO  PLACE  SINGLE  QUOTES  AROUND  THE  FOUR  LETTER ' 
MAP  CODE,  OR  THE  PROGRAM  WILL  NOT  RUN.' 


READ 


AREA 


IF  ( AREA. EQ. 'PAFR')  COAST  =  ' AFR I ' 
IF  (AREA.EQ. 'PASI ' )  COAST  =  'ASIA' 


IF  (ARKA.EQ.'PAUS')  COAST  =  T AUST* 

JF  vAKFA.KQ.  'FKUR' )  COAST  =  'IJURO' 

IF  (AREA.KQ. * PNOR ' )  COAST  =  'NORT' 

IF  (AREA.EQ.'POLJ')  COAST  =  'COAS' 

IF  (AREA.KQ. ' PSOU ' )  COAST  =  'SOUT' 

WRITE  (9,700)  DXMAX ,  DX0R1G,  DXSTP 

700  FORMAT  ( ’ DXMAX  =  ' ,F10.3,5X, 'DXORIG  =  f ,F10.3,5X, 'DXSTP  =  '.F10.3) 
WRITE  (9,701)  DYMAX ,  DY0R1G,  DYSTP 

701  FORMAT  ('DYMAX  =  ' ,F10.3,5X, 'DYORIG  -  ' ,F10.3,5X, 'DYSTP  =  \F10.3) 

C  THE  FOLLOWING  CODE  PRODUCES  A  CONTOUR  PLOT  FROM  THE  DATA  FILE 
C  DOSBRTE  USING  DISSPLA  SUBROUTINES. 

CALL  COMPRS 

CALL  PROJCT('MERCA') 

CALL  YAXANG  (0.) 

CALL  XNAME  ('  ',1) 

CALL  YNAME  ('  ',!) 

CALL  A  REA  21)  (6.,  4.) 

CALL  HEAD  IN  ('I)OSE  RATE  CONTOURS$' ,  100, 1 .5,4) 

ENCODE  (60, 15, LABEL)  Yl.D,  FF 

15  FORMAT  (’YIELD  (KT)  =  \F9.3,'  FISSION  FRACTION  =  ' , F5 . 3 , 

1  1  $ 1 ) 

CALL  IIEAD1N  (LABEL,  100, 1 .  ,4) 

ENCODE  (60, 16, LABEL)  RMED,  BETA l 

1 6  FORMAT  ( ' RMED  ( M I  CRONS )  =  ' , F7 . 3 , '  BETA  =  ' , F6 . 3 , ' $ ' ) 

CALL  HEAD  IN  (LABEL, 100, 1 . ,4) 

ENCODE  (60, 1 7, LABEL)  DLONO,  DLATO 

17  FORMAT  ('GROUND  ZERO:  LONG (DEG. )  =  ',F8.3, 

1  '  LAT(I)EG. )  =  ' ,  F7 . 3 , '  $ '  ) 

CALL  HEAD IN  (LABEL, 100, 1 . ,4) 

CALL  MAPGR ( DXOR I G , DXSTP , DXMAX , DYOR I G , DYSTP , DYMAX ) 

C  CLEAR  OUT  ARRAY  II’KRAY. 

NL  =  LINES'!  ( ll’KRAY, 200,25) 

C  READ  THE  DATA  INTO  THE  PROPER  ARRAYS  AND  CALL  LEGEND-SETTING 
C  ROUTINE. 

CALL  HEIGHT  (0.1) 

I  LINES  =  0 
DO  60  1=1,1000 

READ  (5, ' (3E12.5) ' ,END  =  65)  CX(I),  CY(I),  DRATE(I) 

IF  ( CX ( 1 ) . NE . - 1 . OE+25 )  GOTO  60 

1  LINES  =  I  LINES  +  1 
IF  ( I  LI  NFS. NE. 1 )  THEN 

139 


DSRATE  =  DRATE(L-l) 

WRITE  (9,555)  DSRATE 

555  FORMAT  ('DSRATE  =  *  f FlO-3) 

ENCODE  (25, 20, LEGLAB)  DSRATE 

20  FORMAT  (F10.3,'  ROENTGEN/IIR$ ' ) 

CALL  LLNES  ( LEGLAB, IPKRAY.IL1NES) 

ELSE 

CALI,  LJNES  ('II0TLINE$,,IPKRAY,1LINES) 
WRITE  (9,556) 

556  FORMAT  ( '110TLLNE' ) 

END1F 

60  CONTINUE 
65  CONTINUE 

REWIND  5 


C  GET  THE  WIDTH  AND  HEIGHT  OF  THE  LEGEND  AREA. 

XW  =  XLEGND  ( I PKRAY, 1 LTNES) 

YW  =  YLEGND  ( J PKRAY, ILINES) 

C  SET  UP  A  BLANKED  AREA  FOR  THE  LEGEND. 

YL  =  0.0  -  YW  -  1.0 
XL  =  (6.0  -  XW)  /  2.0 

CALL  BLR EC  (XL,YL,XW  +  0.2, YW  +  0.2,0.02) 

CALL  BLKEY  (ID) 

IF  (AREA.NE. 'USAM' )  CALL  MAPFIL(COAST) 

CALL  MAPI’! L( AREA) 

CALL  TIIKFRM  (.02) 

CALL  FRAME 
CALL  SCLPIC(0. 5') 

WRITE  (9,198)  AREA 
198  FORMAT  ('AREA  =  ' ,A4) 

READ  (5, ' ( 1*10. *5) ' )  YLD 

CALL  RASPLN(0.) 

C  READ  THE  DATA  INTO  THE  PROPER  ARRAYS  AND  CALL  THE  CURVE  PLOTTING 
C  ROUTINE. 

NI’NTS  =  0 
DO  30  1=1,1000 

READ  (5, ' (3E12. 5) ' , END  =  40)  CX(I),  CY(I),  DRATE(l) 

IF  (CX( I ) . EQ.-l .OE+25)  GOTO  35 

NPNTS  =  NPNTS  +  1 

DY(NPNTS)  =  DLA'1‘0  +  CY(I)  /  DEGKM 

DX( NPNTS)  =  DLONO  +  CX(I)  /  (DEGKM*C0S(DY(NPNTS)*P1/180. )) 
WRITE  (9,702)  DX(NPNTS),  DY( NPNTS) ,  NPNTS 
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702  FORMAT  ( ' DX  =  ' .E13.6.5X, 'DY  =  ' ,E13.6,5X, 'NPNTS 

GOTO  30 

35  CALL  CURVE  ( l)X,DY, NPNTS,  1) 

NPNTS  =  0 

30  CONTINUE 

40  CALL  RESET  ('BLNK1') 

CALL  RESET  ( ' BLNK1 ' ) 

CALL  DOT 
CALL  GRID(l,l) 

CALL  BLOFF( LD) 

CALI,  LEGEND  ( IPKRAY,  LI.1NES  ,XL+.  1 ,  YL-h.  J ) 

CALL  ENDPL(O) 

CALL  DONEPL 


100 


CLOSE  (5) 
CLOSE  (9) 
STOP 
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